Unstable eigenvectors and reduced amplitude spaces specifying limit cycles of coupled oscillators with simultaneously diagonalizable matrices: with applications from electric circuits to gene regulation.

Abstract: A fascinating phenomenon is the self-organization of coupled systems to a whole. This phenomenon is studied for a particular class of coupled oscillatory systems exhibiting so-called simultaneously diagonalizable matrices. For three exemplary systems, namely, an electric circuit, a coupled system of oscillatory neurons, and a system of coupled oscillatory gene regulatory pathways, eigenvectors and amplitude equations are derived. It is shown that for all three systems, only the unstable eigenvectors and their amplitudes matter for the dynamics of the systems on their respective limit cycle attractors. A general class of coupled second-order dynamical oscillators is presented in which stable limit cycles emerging via Hopf bifurcations are solely specified by appropriately defined unstable eigenvectors and their amplitudes. While the eigenvectors determine the orientation of limit cycles in state spaces, the amplitudes determine the evolution of states along those limit cycles. In doing so, it is shown that the unstable eigenvectors define reduced amplitude spaces in which the relevant long-term dynamics of the systems under consideration takes place. Several generalizations are discussed. First, if stable and unstable system parts exhibit a slow-fast dynamics, the fast variables may be eliminated and approximative descriptions of the emerging limit cycle dynamics in reduced amplitude spaces may be again obtained. Second, the principle of reduced amplitude spaces holds not only for coupled second-order oscillators, but can be applied to coupled third-order and higher order oscillators. Third, the possibility to apply the approach to multifrequency limit cycle attractors and other types of attractors is discussed.

Introduction

On the one hand, instabilities are at the heart of various self-organization and pattern formation phenomena [1-4]. The way such phenomena enfold in time is typically determined by unstable eigenvectors and eigenfunctions and their corresponding amplitudes. This fundamental framework for the description of systems close to instabilities has been introduced several decades ago [1, 5–8]. Since then, it has produced a steady stream of new applications in contemporary research [9-17] (see in particular recent applications in COVID-19 research [18-20]), which illustrates the important task to make the eigenvector/eigenfunction framework available in research fields that have not yet taken fully advantage of this approach. On the other hand, coupled oscillator systems have been comprehensively studied [21, 22], in particular, to identify their characteristic attractors. A frequently used technique is the phase dynamics approach in which phase dynamics equations are derived from coupled oscillator state equations and stability features of the oscillatory dynamics are primarily discussed on the level of the phase dynamics [22-27]. Alternatively, eigenvalues of coupled oscillator systems are determined to address the nature of the systems’ fixed points and/or extensive numerical simulations are conducted to identify attractors beyond fixed points. The latter approach has been used, for example, to study coupled electric oscillators [28, 29], coupled Fitzhugh–Nagumo oscillators [30-35], and coupled gene regulation oscillators [36-38]. While these approaches provide useful insights, they do not take fully advantage of the aforementioned framework based on unstable eigenvectors and amplitudes | in particular, in the context of systems with simultaneously diagonalizable matrices. Consequently, in the current study, in Sect. 2, we begin with a system of coupled electric oscillators and demonstrate how the concept of unstable eigenvectors and amplitudes can be used to determine the location of and the dynamics on the system’s limit cycle. Subsequently, in Sect. 3, a general framework to address coupled oscillators whose limit cycle attractors are determined by unstable eigenvectors and their amplitudes will be developed. In this context, the theory of simultaneously diagonalizable matrices will play a crucial role. Finally, in Sect. 4, we will work out two further examples from biology: coupled Fitzhugh–Nagumo neurons and coupled gene regulatory oscillators. Sections 2 to 4 will demonstrate that the emerging stable oscillations can be studied with the help of reduced amplitude spaces [1, 13], that is, amplitude spaces that are exclusively spanned by the dominant (unstable) amplitudes. Section 5.1 will generalize this principle for systems exhibiting a coupling between stable and unstable dynamics. It will be shown that under appropriate circumstances the limit cycle dynamics can at least approximately described in terms of an amplitude dynamics the lives in a reduced amplitude space.

Elwakil–Salama coupled electric oscillators

Elwakil and Salama proposed an electric circuit composed of two coupled oscillatory units that can be described by [28]where the pair x, y describes the first unit and u, v describe the second unit; x and u are voltage-like variables, whereas y and v are current-like variables. However, in Eqs. (1)–(4), re-scaled variables are shown that correspond to dimensionless quantities. The parameters and correspond to the quality factor of the electric circuit (Q) and the gain (G) of the nonlinearity that connects the two electric units. Q and G, again, are dimensionless quantities. The angular frequency parameter can be eliminated by rescaling time. However, in the context of numerical simulations, it is helpful to keep in the equations. Equations (1)–(4) can be equivalently expressed as the coupled system of second-order equationswith . The transformation into Eq. (5) is a key step, because it allows us to identify two new variables, one of which is the aforementioned variable p. The second new variable is given by . From Eq. (5), it follows that the evolution of s is determined by:and the eigenvalueswith . In Eq. (7) and in what follows i is the imaginary unit. Equation (6) describes a damped oscillation. In contrast, for from Eq. (5), it follows that:Equation (8) describes a van-der-Pol oscillator [39] of the formwith a potential function , and . Linearizing Eq. (8) at the fixed point yieldsFrom Eq. (10), it follows that the evolution of p close to the fixed point is characterized by the eigenvalues:with , . The eigenvalues shown in Eq. (7) and (11) have been previously derived using a different approach [28]. Let us evaluate Eq. (11). For , we have and . For , we have and . For , we have and withand . For , we have , such that becomes real with . For , the eigenvalues remain real valued.

In the context of the current study, the parameter domain is of particular interest. At with , a Hopf bifurcation takes place and p describes a self-oscillator. This also implies that the original coupled oscillator system (1)–(4) exhibits a Hopf bifurcation at . From the definitions of s and p, it follows that the 2D state vector can be expressed like:withAs will be shown in Sect. 3, the vectors and correspond to non-normalized eigenvectors of appropriately defined matrices of the electric circuit defined by Eqs. (1)–(4). For , the vector corresponds to the stable eigenvector (because of ), whereas corresponds to the unstable eigenvector (because of ). Equation (13) expresses the state vector in terms of a superposition involving the new coordinates s and p. In this context, s and p are frequently referred to as amplitudes [1, 13, 18]. From Eqs. (6) and (13), it then follows that:Equation (15) states that determines an axis that specifies in the x–u state space the orientation of the stable limit cycle of the electric circuit. That is, the relevant amplitude space is reduced from 2D (s and p) to 1D (p only).

Let us illustrate the considerations made on the Elwakil–Salama circuit by means of a numerical simulation. We solved Eqs. (1) to (4) for , , s and initial conditions , , , and . The differential equations were solved by means of a standard Euler forward scheme with fixed time step of 0.001s. Figure 1 shows x,u (panel (a)) and s, p (panel (b)) as functions of time. As expected, the dynamics evolved towards a limit cycle with and . Importantly, Fig. 2 demonstrates the role of the unstable eigenvector . Panels (a) and (b) show the phase curve u(x) in the x–u state space. Panel (a) shows u(x) for the interval to s. Panel (b) shows the full-phase curve. is shown as well. As can be seen, the phase curve converged towards the direction of when time elapsed. As such, the circuit defined by Eqs. (1)–(4) can be described in the 4D space x, y, u, v or alternatively in the 4D space . Panel (c) of Fig. 2 shows the simulated trajectory in the 3D subspace . The direction specified by the unstable eigenvector is indicated in a x–u plane, as well. In panel (c), the stable limit cycle can be clearly seen. As expected, we found that the limit cycle was oriented in the 3D space along . In other words, the limit cycle was located in a plane determined by .

Fig. 1

Simulation results (part 1) obtained for the Elwakil-Salama circuit (1)–(4). Panel a: x (solid line) and u (dotted line) as functions of time. Panel b: s and p over time. The figure demonstrates for the Elwakil–Salama circuit that while in state space, two variables (x and u) are needed to characterize the limit cycle dynamics of the electronic system, in amplitude space, only a single variable (p) is needed. For parameters and initial conditions, see text

Fig. 2

Simulation results (part 2) obtained for the Elwakil–Salama circuit (1)–(4). Panels a and b show simulated phase curves u(x) (solid black lines) and the direction (dashed gray lines) defined by the unstable eigenvector in the 2D x–u state space. Panel a shows a detail of panel b for the first 19 seconds of the simulation. Panel c shows the simulated trajectory (solid black line) of the electric circuit in the 3D space spanned by and the direction (dashed gray line) of illustrated in a x–u plane. Panels a–c: The squares indicate the initial state used in the simulation. The figure demonstrates that it is the unstable eigenvector (or order parameter) that characterizes the orientation of the limit cycle. The figure illustrates graphically Eq. (15)

General coupled oscillators with limit cycles specified by unstable eigenvectors

In this section, a general framework of coupled oscillators with limit cycles specified by unstable eigenvectors is developed. To this end, we first return to the example of the Elwakil–Salama circuit. As shown in Sect. 2, the circuit (1)–(4) can be expressed in the 4D state space x, u, by means of Eq. (5). Linearizing Eq. (5) at , , the linearized evolution equations can be written likeImportantly, Eq. (16) involves two matrices. Both matrices exhibit the same eigenvectors, namely, and . They correspond to the non-normalized eigenvectors reported in Eq. (14). Since the Elwakil and Salama model exhibits as unstable eigenvector , from Eq. (15), it follows that the stable limit cycle of the Elwakil–Salama circuit is determined by the ”anti-diagonal axis” in the x–u space (see also Fig. 2).

Let us derive a general framework of coupled oscillator systems for which the unstable eigenvector can point in an arbitrary direction in the x–u state space, such that those systems can exhibit stable limit cycles that are oriented in arbitrary directions. The following considerations apply to physical, chemical, and biological systems in general and are not restricted to electric circuits. Let x(t) and u(t) denote the state variables of two coupled units that satisfyand exhibit the fixed point , . Let denote the 2D state vector. Then, the linearized version of Eq. (17) at the fixed point can be written likewhere A and B are matrices. Next, we require that the matrices A, B are simultaneously diagonalizable. Simultaneously, diagonalizable is the property when two matrices are diagonalizable and exhibit the same eigenvectors [40]. That is, A, B exhibit the property that there exists a diagonalization matrix M, such thatwhere and are diagonal matrices with diagonal elements , and , , respectively. Moreover, let and denote the column vectors that constitute M likeSince (by assumption) exists, the two vectors and are linearly independent. Then, it follows that and are the eigenvectors of A and B and are the corresponding eigenvalues with and like (e.g., ). With the help of Eqs. (19) and (20), Eq. (18) can be transformed into diagonal form (see appendix 7.1)The eigenvalues of the two linear dynamical systems (s and p) defined in Eq. (21) are given byFinally, the state can be expressed in the basis of the eigenvectors and with the help of the variables s and p likeAs indicated in Sect. 2, we will refer to s and p as amplitudes and the 2D s–p space as amplitude space [1, 13, 18]. Equation (23) describes the mapping from amplitude to state space.

In what follows, let us consider a Hopf bifurcation scenario. Let us assume that there exists a bifurcation parameter in an interval with a critical value . Let us assume that the coupled system exhibits the property that holds for any , whereas , , and . Finally, let us further assume that holds for . Then, a Hopf bifurcation takes place at in the subsystem described by the variable p. Moreover, the dynamics of the amplitude s exhibits a stable fixed point. Therefore, the Hopf bifurcation in the p dynamics carries over to the original 4D system defined by Eq. (17). For , Eq. (17) exhibits a Hopf bifurcation. The Hopf bifurcation is characterized by the two eigenvectors and . The vector describes a stable eigenvector (because of ), whereas describes an unstable eigenvector for (because of for ). For appropriate nonlinear terms in Eq. (17), we re-obtain Eq. (15), which is here repeated asFurthermore, if for appropriate nonlinear terms, the p-dynamics describes a stable limit cycle dynamics, then the 4D system exhibits a stable limit cycle whose orientation is determined by and dynamics is determined by the amplitude p. The relevant amplitude dynamics takes place in a 1D space, that is, in a reduced amplitude space.

While the linear terms of the model allow us to identify the eigenvalues and eigenvectors and the stability of the fixed point, the nonlinear terms determine the approach towards the limit cycle attractor. Therefore, let us turn next to the nonlinear terms. To this end, let us re-write Eq. (17) in the formwhere contains only nonlinear terms. Let us consider the special case in which s(t) satisfies the linear dynamics described in Eq. (21). That is, we need to add nonlinear terms only to the evolution equation of p shown in Eq. (21). This implies that shows in the direction of likeIn the simplest, non-trivial case depends on p and only, such thatFor example, let us consider the van-der-Pol term and the Rayleigh term [39] . In this case, from Eq. (27), we obtain as an intermediate resultTo obtain an explicit, final form of the coupled oscillator system, we need to express p in terms of x and u using Eq. (96). To this end, we may use the bi-orthogonal vectors and [1, 18, 41] (for details see appendix 7.2).

From our considerations above, it follows that the dynamics of s and p is given by:Assuming that the assumptions and at of the aforementioned Hopf bifurcation scenario hold, then Eq. (99) describes a coupled oscillator system that exhibits a stable limit cycle determined by Eq. (24) that is characterized by the unstable vector , on the one hand, and the van-der-Pol-Rayleigh dynamics of p(t), on the other hand. The latter is defined by Eq. (30).

In general, and are (by definition) linearly independent vectors (see above). In the special case when and constitute a pair of orthogonal vectors, the matrix M assumes the form of the 2D rotation matrixinvolving the angle , such thatSince the basis is an orthogonal basis, the bi-orthogonal vectors are given by and .

First, let us illustrate by means of a numerical simulation the special case of an orthogonal basis. We solved numerically (Euler forward scheme with fixed single time step of 0.001s) the coupled oscillator model (99) for /s, /s, /s, /s, , and s. Consequently, the eigenvalues were per second and per second, such that and corresponded to the stable and unstable eigenvectors, respectively. Explicitly, the vectors with were given byThe vectors (33) were used as coefficients of the matrix M (see Eq. (31)). The matrices A and B occurring in Eq. (99) were then computed from the diagonal matrices and (whose coefficients are listed above) by inverting the steps shown in Eq. (19). That is, and . The bi-orthogonal vector in Eq. (99) was simply given by as discussed above. The following initial conditions were used: , , , and . Figures 3 and 4 present the simulation results. Panel (a) of Fig. 3 shows x and u as functions of time. It can be seen that the dynamics converged towards a stable limit cycle dynamics. Panel (b) shows the variables s and p over time. As expected, s converged towards zero, whereas p approached a stable limit cycle dynamics. Figure 4 shows the oscillator system dynamics in 2D (panels (a) and (b)) and 3D (panel (c)) state spaces. Panel (a) show the phase curve u(x) for the first 11 s, whereas panel (b) shows the entire simulated trajectory. The directions specified by vectors and are shown in panels (a) and (b), as well. By construction, and define axes at angles of 100 and 10, respectively, with respect to the horizontal axis (i.e., the x axis). As can be seen in panels (a) and (b), after a transient period, the phase curve u(x) approached the axis defined by the unstable eigenvector . Subsequently, the dynamics took place along the axis specified by as predicted by Eq. (24). Panel (c) reveals the shape of the limit cycle in the 3D space given by x,u and . is drawn in an x–u plane, as well. Panel (c) demonstrates that in the 3D space spanned by x, u, and , the limit cycle is located in a plane specified by the axis. Note that (as mentioned above), we used the simulation parameters and . Consequently, the limit cycle dynamics was determined by a Rayleigh limit cycle oscillator.

Fig. 3

Simulation results (part 1) of the coupled oscillator model (99) in the case of the orthogonal basis given by Eq. (33). Panel a: x (solid line) and u (dotted line) as functions of time. The horizontal bar starting at s reflects , where is the maximal value of x obtained numerically during the simulation period and (see also Eq. 34). Panel b: s and p over time. The figure demonstrates for a generic system of two coupled nonlinear oscillators (17) with state variables x and u that under appropriate conditions, only a single variable, namely, the unstable amplitude p, is needed to characterize the limit cycle dynamics. For parameters and initial conditions, see text

Fig. 4

Simulation results (part 2) of the coupled oscillator model (99) in the case of the orthogonal basis given by Eq. (33). Panels a and b show simulated phase curves u(x) (solid black lines) and the axes given by the eigenvectors and (dashed gray lines) in the 2D x–u state space. Panel a shows a detail of panel b for the first 11 seconds of the simulation. Panel c shows the simulated trajectory (solid black line) of the coupled oscillator system in the 3D space spanned by . The axis of (dashed gray line) is plotted in a x–u plane as well. Panels a–c: The squares indicate the initial state used in the simulation. The figure demonstrates that the orientation of the limit cycle is entirely determined by the unstable eigenvector and, in doing so, exemplified and visualizes Eq. (24)

From Eq. (24), it follows that on the limit cycle:holds with . From Eq. (33) for the simulation shown in Figs. 3 and 4, we get . That is, the oscillation amplitude of u should be about 5 times smaller as compared to the amplitude of x. In panel (a) of Fig. 3 the maximal value of u predicted from is shown as horizontal bar (starting at s). We found that the numerical maximal values of u(t) fitted very well the theoretical ones.

Next, let us illustrate the general case of a non-orthogonal basis by means of a numerical simulation (to this end, again, an Euler forward scheme with time step 0.001s was used). In this case, the model involves bi-orthogonal vectors as discussed in the appendix 7.2. We solved Eq. (99) (see again appendix 7.2) for the same diagonal matrices and as in the previous example. That is, we used /s, /s, /s, and /s. Consequently, the eigenvalue were again given by per second and per second. We also considered again the Rayleigh oscillator case with and s. However, in this second simulation, the eigenvectors were not defined by Eq. (32). Rather, for the simulation, we usedThe vectors are linearly independent, but they are not orthogonal. The same initial conditions were used as in the previous simulation. From Eqs. (29) and (30), it then follows that the s and p amplitude dynamics of this second simulation should be identical to the s and p amplitude dynamics of the previous simulation. In contrast, the dynamics of the state variables x and u should differ across the two simulations. The matrices A and B of Eq. (99) were computed again from and with M defined by Eq. (20). Moreover, of Eq. (99) was taken as the second row vector of as discussed above (or see again Refs. [18, 19]).

Figures 5 and 6 present the simulation results. The panels of Figs. 5 and 6 present the same kind of variables and quantities as the panels of Figs. 3 and 4. As can be seen in Fig. 5, as expected, x(t) and u(t) approached a limit cycle dynamics with . In particular, as discussed above, we found that the dynamics of s and p in this second simulation was identical to the s and p dynamics of the previous simulation [compare panels (b) of Figs. 3 and 5]. Moreover, for this second simulation, Eq. (34) predicts that on the limit cycle, the peak values are 4 times larger than the peak values of because . We found that this was indeed the case [see the horizontal bar in panel (a) of Fig. 5]. Moreover, u and x exhibited an anti-phase synchronization consistent with . Panel (a) of Fig. 6 illustrates the non-orthogonal basis spanned by and . As predicted by Eq. (24), irrespective of the fact that and were non-orthogonal, the simulated trajectory converged over time towards the axis defined by . Panel (c) depicts the shape of the limit cycle in the 3D space given by x, u and . For comparison purposes, is drawn in an x–u plane of the 3D space. Just as for the previous simulation, panel (c) of Fig. 6 demonstrates for the second simulation that the limit cycle in the 3D space spanned by x, u, and is located in a plane specified by the axis.

Fig. 5

Simulation results (part 1) of the coupled oscillator model (99) in the case of the non-orthogonal basis given by Eq. (35). Panel a: x (solid line) and u (dotted line) as functions of time. The horizontal bar starting at s reflects and predicts . Panel b: s and p over time. For parameters and initial conditions, see text

Fig. 6

Simulation results (part 2) of the coupled oscillator model (99) in the case of the non-orthogonal basis given by Eq. (35). Panels a–c of Fig. 6 show the same kind of variables and quantities as the respective panels in Fig. 4. As far as panel a is concerned, u(x) for the first 5 s of the 50 s simulation is shown

Our approach applies to nonlinear systems whose linearized evolution equations exhibit simultaneously diagonalizable matrices. Let us briefly dwell on implications of this requirement. Obviously, on the one hand, this constraint implies that the approach does not apply to all kind of dynamical systems. However, on the other hand, it does not necessarily require that we discuss coupled identical systems as it was the case for the Elwakil–Salama model (see Eqs. (1)–(4)). Let us illustrate this point for coupled second-order systems that satisfy Eq. (17). The corresponding linearized evolution equation in the most general case involves 8 independent coefficients and with that correspond to the matrix elements of A and B, respectively. If we deal with two coupled identical oscillators, then Eq. (18) reads explicitlyand involves only 4 independent coefficients. In contrast, if the two matrices A and B occurring in Eq. (18) can be simultaneously diagonalized, then they can exhibit up to 6 independent coefficients. These coefficients can be expressed in terms of the 4 diagonal elements , , , , the independent vector element and the independent vector element . Note that and likewise are not independent coefficients due to the normalization constraint. In other words, from Eq. (19), it follows that A and B can expressed likeIn conclusion, the systems that can be addressed by our suggested approach are more general than systems composed of coupled identical units but less general than the most general case of systems. Let us exemplify these considerations for the numerical simulations presented above. For the first model whose solutions are displayed in Fig. 3 and 4 and involves the orthogonal basis (33), the linearized model (18) reads explicitlyThe two coupled units exhibit the same off-diagonal coefficients (i.e., the same coupling coefficients). However, they exhibit different diagonal elements. For the second model whose solutions are displayed in Fig. 5 and 6 and involves the non-orthogonal basis (35), the linearized model (18) readsAll eight coefficients differ from each other, which demonstrates again that the approach does not require that we consider coupled identical units. Note that while the numerical values of all eight coefficient differ, the maximal number of independent coefficients of the linearized model is six. That is, the eight values shown in Eq. (39) have been computed from the six independent values , , , , , and listed above.

Further applications

Coupled Fitzhugh–Nagumo neurons

The Fitzhugh–Nagumo model is a limit cycle oscillator model that describes periodically firing neurons [39, 42, 43]. That is, each round-trip through the limit cycle is interpreted as a spiking of a neuron and a subsequent refractory phase of the neuron. Coupled Fitzhugh–Nagumo neurons have been discussed in the literature in the context of synchronized and chaotic brain activity (see, e.g., Refs. [30-33]) and, in particular, in the context of chimera states [34, 35]. A linear coupling between the two neurons has frequently used. In what follows, two coupled Fitzhugh–Nagumo neurons are considered, where x, y denote the state variables of the first neuron and u, v denote the variables of the second neuron. The coupled neuron model readsand involves a linear coupling term with . The function F is frequently taken as a cubic nonlinearity of the form with [39]. In this context, a may be consider as some kind of pumping parameter. In contrast, c may be interpreted as damping parameter (see also the eigenvalue analysis below). As far as the measurement units of the model variables are concerned, the variables x, u may be regarded as voltage-like variables measured in Volt, while y, v may be considered as current-like variables measured in Ampere. However, frequently, Fitzhugh–Nagumo neurons are model by means of dimensionless variables. Accordingly, in what follows, the state variables x, y, u, v will be considered as dimensionless quantities. The model defined by Eqs. (40)–(43) can be equivalently expressed as Equation (44) is an important intermediate result, because it shows that the coupled neuron model assumes the general form (25) and it allows us to introduce two new variables, which we will do next. Introducing the variables and , from Eq. (44), we obtainwhere and are nonlinear terms. From Eq. (21) and (22) and the linear parts in Eq. (45), it follows that the eigenvalues read:For the current study, it is sufficient to consider the case discussed by Jenkins [39] with and . The latter condition implies that the expression occurring in Eqs. (45) and (46) is positive. In this case, from Eq. (46), it follows that forthe eigenvalues are either real and negative or they are complex-valued and exhibit negative real parts. Moreover, are either real and positive or they are complex-valued and exhibit positive real parts. Consequently, the variable s describes a stable dynamics towards (at least in the linear domain), whereas the variable p describes an unstable dynamics. The inequality (47) refers to a coupled neuron system for which the pumping parameter a is sufficiently large relative to the damping parameter c, such that holds. However, the coupling between the neurons as measured by D is sufficiently strong relative to the individual pumping a of the neurons, such that holds. In this case, from , it follows that for appropriately chosen model parameters, the variable s converges to zero. If so, the dynamics of p(t) is determined bywhich is the equation of a van-der-Pol oscillator with potential . The 2D state vector can be expressed likewithThe vectors correspond to the eigenvectors of the simultaneously diagonalizable matrices A and B of the linear parts of Eq. (44), and in view of Eq. (49), s and p are considered as the amplitudes of the coupled neuron model. Since for holds under the aforementioned conditions, Eq. (49) leads to Eq. (24) again. Accordingly to Eq. (24), the orientation of the stable limit cycle is determined by the unstable (non-normalized) eigenvector as defined in Eq. (50). The dynamics on the limit cycle is determined by Eq. (48).

To illustrate these properties of two coupled Fitzhugh–Nagumo neurons, we simulated the model (40)–(43) (with the help of an Euler forward scheme with time step 0.001s) for the parameters /s, /s, and /s, such that Eq. (47) was satisfied with and the eigenvalues were given by per second and per second. Furthermore, we used the initial conditions , , , and . Figures 7 and 8 present the simulation results. Figure 7 presents x,u (panel (a)) and s, p (panel (b)) as functions of time. The neuron dynamics evolved relatively quickly towards a limit cycle with and . Panels (a) and (b) of Fig. 8 present the stable and unstable axes defined by and and the neural dynamics in terms of the phase curve x(u) for the first 7.5 s (panel (a)) and the entire simulation duration (panel (b)). As can be seen, the simulated neural dynamics approached the axis given by , which is consistent with Eq. (24). Panel (c) of Fig. 8 presents the simulated trajectory in the 3D subspace . The stable limit cycle can be clearly seen. The limit cycle was oriented along the direction specified by .

Fig. 7

Simulation results (part 1) of the coupled neuron model (40)–(43). Panel a: x (solid black line) and u (dotted gray line) as functions of time. Panel b: s and p over time. The figure illustrates that while from a state space perspective (panel a), the model describes synchronized neural activity in terms of two synchronized neural signals x and u, the amplitude space perspective (panel b) uncovers that in fact the coupled system exhibits a single dominant variable p (the unstable amplitude) that ”drives” and determines the synchronization of x and u observed in the state space. For parameters and initial conditions, see text

Fig. 8

Simulation results (part 2) of the coupled neuron model (40)–(43). Panels a and b show simulated phase curves u(x) (solid black lines) and the axes given by the eigenvectors and (dashed gray lines) in the 2D x–u state space. Panel a shows a detail of panel b for the first 7.5 s of the simulation. Panel c shows the simulated trajectory (solid black line) of the coupled neuron model in the 3D space spanned by . The axis of (dashed gray line) is plotted in a x–u plane as well. Panels a–c: The squares indicate the initial state used in the simulation. The figure demonstrates that the neural two units’ system exhibits an underlying unstable amplitude (or order parameter) that determines the orientation of the limit cycle in state space. In doing so, the figure exemplifies Eq. (24) for the concrete application of two coupled Fitzhugh–Nagumo neurons (40)–(43)

Coupled Goodwin oscillators: synchronized inter-cellular gene regulation

One of the earliest models that describes oscillatory gene regulation within a single cell is the Goodwin oscillator [44-46]. The three-variable Goodwin oscillator reads [44, 45]with and . In Eqs. (51)–(53), the variables , , denote re-scaled quantities and describe a gene regulation loop in terms of the amount of mRNA () produced by a DNA promoter, the amount of inactive proteins () produced from the mRNA, and the amount of active (or activated) proteins () obtained from the inactive proteins. The active proteins regulate the promoter via the function f. In what follows, a Hill-type function f of the form: [45, 47, 48]with and will be used, where n is an integer. It is possible to rescale the variables even further, such that and [45]. However, for the purpose of conducting numerical simulations, it is more convenient to keep the parameters in the analytical considerations. As such, k can be used to adjust the time scale of the dynamics. The parameters B and determine the fixed point of the regulatory unit (see below). Following earlier work in which the Goodwin model has been studied in detail (see e.g., Refs. [38, 45, 47]), we will briefly analyze some key features of the particular version of the Goodwin model defined by Eqs. (51)–(54). The fixed point of the model (51)–(54) is given byEquations (51)–(53) can be written as a third-order differential equation of the formLet denote deviations from the fixed point value like . Then, the linearized version of Eq. (56) readswithUsing Eq. (55), we obtainThe eigenvalues of Eq. (57) are given by andConsequently, for , , , the real parts of are negative, zero, and positive, respectively. This implies that a Hopf bifurcation takes place at with . From Eq. (55), it follows that X increases monotonically as a function of . Likewise, from Eq. (59), it follows that G increases monotonically as a function of X. Consequently, G increases monotonically as function of the ratio and the ratio may be considered as bifurcation parameter. Note that from Eq. (59) and the requirement for an unstable fixed point, it follows that . Since holds, the inequality cannot be satisfied for . The smallest possible integer that allows for is , which is a key observation that has been made in the previous literature [38, 45, 47].

Coupled Goodwin oscillators have been considered in various studies. In his benchmark study, Goodwin [44] considered the gene regulation within cells by means of two coupled pathways. More recent, gene regulation across cells (i.e., inter-cellular interactions between cells) have been modeled by means of coupled Goodwin oscillators [36-38]. In line with those studies, in what follows, two cells are considered that produce proteins (cell 1) and (cell 2), respectively. The active forms of the proteins are given in the amount (cell 1) and (cell 2). Previous work suggest that when the protein produced by one cell act as repressors of the promoter activity of the other cell, then the two cells under appropriate conditions can synchronize their oscillatory activity [37]. In what follows, the approach presented in the previous sections will be used to examine this proposed mechanism of synchronized inter-cellular gene regulation. The aim is to illustrate the basic idea. For this purpose, it is sufficient to consider a simplified model in which the uncoupled gene regulation loops of the proteins and exhibit the same fixed points . Moreover, the repressor impact h from one cell onto the other cell will be described by a linearized function and is assumed to leave the fixed point unaffected. This implies that the feedback function that describes the interactions between the cells readswith . Accordingly, the gene regulatory system composed of the two coupled cells reads The following two steps are conducted to introduce two new variables s and p, again, as in the previous examples. First, the model is written in terms of third-order differential equations that read Second, we introduce the deviations and defined by and , respectively. Then, the coupled linearized equations readLet us define s and p like and . Then, the linearized evolution equations of s and p are given byThe eigenvalues read (assuming ), andIn the context of the current work, the case of interest is given byIn this case, the eigenvalues satisfy and , such that the variable s describes a stable dynamics approaching (at least as long as the linearized equations hold), whereas the variable p describes an unstable dynamics characterized by an initial exponential increase of p. The state vector can be expressed likewithand motivates to consider s and p as amplitudes. If (in view of the negative real parts of the eigenvalues and ) for holds, then, in analogy to Eq. (24), from Eq. (73), it follows that:Accordingly, the orientation of the stable limit cycle is determined by the unstable (non-normalized) eigenvector . Substituting into the nonlinear evolution equation of p(t) [that can be derived from Eq. (68)], we conclude that in this case, the dynamics on the limit cycle is determined by

To illustrate these considerations, we solved Eqs. (62)–(67) for , , , and which implies and such that the inequality (72) was satisfied like . Since the Goodwin oscillator has been frequently discussed in the context of circadian rhythmic activity on the cellular level that is autonomous (and not driven by the day–night cycle) [36, 46, 49], we selected the parameter k as /d, such that the model produced oscillations with a period of approximately 1 day. Finally, the following initial conditions were used: and , , . For simulation purposes, an Euler forward scheme was used with fixed time step of 0.001d. Figures 9 and 10 present the simulation results.

Fig. 9

Simulation results (part 1) of the coupled Goodwin oscillator model (62)–(67). Panel a: (solid black line) and (dotted gray line) as functions of time. Panel b: s and p over time. The figure illustrates for the model (62)–(67) the existence of a single dominating amplitude p (panel b) that underlies the synchronized limit cycle gene regulatory activity observed on the state variable level (panel a). For parameters and initial conditions, see text

Fig. 10

Simulation results (part 2) of the coupled Goodwin oscillator model (62)–(67). Panels a and b show simulated phase curves (solid black lines) in the 2D – state space and the unstable eigenvector axis (dashed gray lines). Panel a shows a detail of panel b for the first 5 days of the simulation. Panel c shows the simulated trajectory (solid black line) of the model in the 3D space spanned by . The axis of (dashed gray line) is plotted in a – plane, as well. Panels a–c: The squares indicate the initial state used in the simulation. The figure demonstrates that irrespective the fact that the coupled Goodwin oscillator model involves third-order oscillators and not second-order oscillators as in the previous examples, the model exhibits an unstable eigenvector () that specifies a direction (panels b) and a plane (panel (c)), in which the gene regulatory limit cycle activity can be found

Figure 9 presents the solutions , (panel (a)) and s(t), p(t) (panel (b)). As expected, the gene regulatory system evolved towards a limit cycle reflecting a synchronized state with and . Panels (a) and (b) of Fig. 10 present the unstable axis given by and the phase curves for the first 5 days (panel (a)) and the entire simulation period of 10 days (panel (b)). We found that the simulated dynamics approached the axis given by as predicted by Eq. (75). Panel (c) of Fig. 10 presents the simulated trajectory in the 3D subspace together with the direction specified by . We found that the stable limit cycle was clearly oriented along .

While Fig. 10 demonstrates the spatial aspect of Eq. (75) about the orientation of the limit cycle in the x–u space, let us illustrate more explicitly the dynamic aspect of Eq. (75). Accordingly, the dynamics in the reduced amplitude space given by the unstable amplitude p entirely determines the dynamics on the limit cycle. To this end, let us compare p(t) as obtained from with p(t) as computed from Eq. (76). Equation (76) was solved numerically in terms of three coupled first-order differential equations with , , and , which are related to the original state variables like , , and . We solve the original model (62)–(67) for the same parameters and initial conditions as used in the previous simulation but for a longer duration of s to make sure that the dynamics approached the stable limit cycle. Furthermore, we solve Eq. (76) in terms of , , , for the same parameters and initial conditions and the same simulation duration. Figure 11 shows the stable limit cycle in the space –p space (i.e., – space) as obtained from the last 10s interval of the 30 s simulation of (62)–(67) and as obtained from the same last 10 s interval of the simulation of Eq. (76). As can be seen in Fig. 11, the solutions obtained by the two methods were identical, as expected.

Fig. 11

Comparison of the limit cycle dynamics obtained from the full Goodwin model (62)–(67) and the amplitude equation (76): versus p(x(t), u(t)) computed from Eq. (62)–(67) (solid thick black line) and versus p(t) computed from Eq. (76) (dotted gray line) are shown. The first 20s of the 30 s-long simulation were removed. Only the dynamics during the last 10 s interval is shown. See text for details. The figure illustrates that the closed evolution equation for the unstable amplitude (76) is sufficient to describe the evolution of the system along its limit cycle. In doing so, the figure also provides a graphical illustration of Eq. (75)

Generalizations

When s does not vanish but corresponds to a “fast” variable

In the previous sections, systems have been discussed for which the unstable eigenvector and the associated amplitude p entirely determine the systems’ limit cycles. In the case of more general systems, both stable and unstable eigenvectors and their amplitudes contribute to the overall emerging pattern or organization. Typically, the unstable eigenvectors and their amplitudes make the major contributions, while the stable eigenvectors and their amplitudes make minor contributions [1, 13]. Nevertheless, these minor contributions do not vanish. Of particular interest is the case when the stable amplitudes evolve relatively quickly as compared to the unstable amplitudes [1, 13]. The following considerations exemplify this situation for the general coupled oscillator model outlined in Sect. 3. While in Sect. 3, examples have been addressed for which the dynamics of the stable amplitude s does not depend on the unstable amplitude p, we consider next the case in which the dynamics of s is coupled to p. Let us consider the fundamental example of a negative-feedback loop in which the unstable amplitude ”pumps” the stable amplitude, while the stable amplitude leads to a damping of the unstable amplitude [13, Chap 4.]. Accordingly, in Eq. (29), we add on the right-hand-side a driving force p, while in Eq. (30), we replace the van-der-Pol damping term with a mixed term , such that the amplitude equations readwith , , , and . In terms of the state variables , we consider a model that assumes the form (25). To simplify the presentation, we consider the special case of a system involving orthogonal eigenvectors (the general case that involves non-orthogonal vectors can be worked out in analogy to the example given in appendix 7.2). In this case, the expressions p and can be computed from and . Consequently, on the level of the state variables x, u, , and , the amplitude equation model given by Eqs. (77) and (78) readsAs such, the solution of Eq. (77) can be found with the help of a Green’s function (or response function) G(z) [24, 50] like and , where and denote the homogeneous and particular solutions, respectively. However, s is assumed to evolve on a short time scale relative to p, which implies that decays to zero relatively quickly and can be neglected. As far as is concerned, since s describes a damped harmonic oscillator driven by p (see Eq. (77)), s follows p with a short delay . When assuming that oscillations of p are relatively small in magnitude, such that p(t) describes an approximately sinusoidal oscillation, then s can be expressed in terms of p likewhere f is proportional factor. Note that in Eq. (80), the small delay approximation has been used [51]. It is well known that the proportionality factor f and delay are given by [52]andwith . Moreover, the parameter occurring in Eqs. (81) and (82) corresponds to the angular driver frequency, that is, the angular frequency of the variable p assuming the sinusoidal approximation holds. In this case, . The variable s can be considered as fast variable (relative to p) if and (i.e., ) holds. In particular, implies that vanishes relatively quickly, while implies that (see Eq. (82)). Similar to the adiabatic elimination technique and other techniques used to study fast variables in other contexts [1, 13, 53, 54], based on the previous considerations, an approximative amplitude equation model can be constructed. Let and denote the approximations of s and a, respectively (i.e., and ). Then, is obtained by a mapping from and like (see Eq. (80))and satisfies the closed evolution equation (see Eq. (78) with s replaced by )Equation (84) provides an approximative evolution equation for the amplitude dynamics of the oscillator model. Let us turn next to an approximative description of the orientation of the limit cycle. As mentioned at the beginning of this section, the theory of pattern formation states that an emerging state given in terms of a superposition of eigenvectors may be composed both of stable and unstable eigenvectors, where the unstable ones typically make the major contributions [1, 13]. Consequently, when s does not vanish, then on the limit cycle, the state as defined by Eq. (23) will not be specified entirely by (i.e., Eq. (24) will not hold). If is sufficiently small and p(t) evolves relatively slow (as assumed), such that in a first-order approximation can be neglected versus p, then from Eq. (80), we obtain . Substituting this approximation into Eq. (23), we obtainThat is, after an initial period during which becomes negligibly small, the 2D state vector satisfieswhere denotes an effective vector that specifies the approximative orientation of the limit cycle in the x–u state space.

Fig. 14

As in Fig. 13 but for the ”late” period from 80 to 100s

Let us illustrate these considerations by a numerical simulation. Just as in Sect. 3, we considered the state variables as dimensionless quantities. We solved numerically the coupled oscillator model (79) for /s, /s, /s, /s, /s and /s. We used an orthogonal basis (33) with . Consequently, the model exhibited the eigenvalues per second and per second with five times larger in the amount than . Moreover, /s and /s, such that was three times larger than . Note that from , it follows that the oscillator system had a characteristic period of s. The matrices A and B were computed as in Sect. 3. The initial conditions , , , and were used. The numerical simulation was conducted by means of an Euler forward scheme with fixed time step of 0.001s.

Fig. 13

Amplitude dynamics of the coupled oscillator model (79). Panel a: s(x, u) and p(x, u) computed from Eq. (79) are shown as thick black lines and s(t) and p(t) computed from Eq. (77) and (78) are presented as thin gray lines. Since the solutions were identical the thin gray lines can be seen on top of the thick black lines. Panel b: s(x, u) and p(x, u) computed from Eq. (79) (as in panel a) are shown again as thick black lines. and computed from Eq. (83) and (84) are presented as thin gray lines. The gray lines deviated slightly from the black lines, which can be seen more systematically in Fig. 14. The figure demonstrates that the closed approximative evolution equation (84) for the unstable amplitude yields a good approximation to the exact solution. Panels a and b: all solutions are shown for the initial period of 80s

Figures 12, 13, 14 and 15 present the simulation results. Figure 12 shows x(t) and u(t) for the initial period of 80s (panel (a)) and for the ”late” period from 80 to 100s (panel (b)). In the ”late” period, the system dynamics evolved approximately on the limit cycle. Figures 13 and 14 demonstrate the use and goodness of the approximation via Eqs. (83) and (84). Figure 13 shows the solutions s(t) and p(t) as computed from the state variables x and u via and , respectively. These solutions are considered as the exact solutions and are presented both in panels (a) and (b). Panel (a) of Fig. 13 shows the solutions s(t) and p(t) of the exact amplitude equations (77) and (78) (that we solved numerically) as well. As expected, the solutions were identical to those obtained from the state variables. In contrast, panel (b) of Fig. 13 shows the solutions and of the approximative amplitude equations (83) and (84) (that were solved again numerically). We found that the functions and were a good fit to s(t) and p(t) but deviated slightly from the exact solutions. Figure 13 shows these results for the initial period of 80s. Figure 14 present the same analysis results for the ”late” period that reflects the evolution on the limit cycle. In particular, panel (b) reveals that and deviated slightly from the exact solutions at the extreme values that can be seen at s, s, and s. Figure 15 shows the system’s dynamics as phase curves in the 2D x–u space (panel (a)) and the 3D space given by x, u and (panel (b)). The entire 120s dynamics is shown in both panels. Panel (a) presents the axes defined by , , and . The axis forms a angle with the horizontal axis (i.e., the x axis). As such and are orthogonal to each other. However, this property cannot be seen in panel (a), because the x and u axes do not show the same ranges. As can be seen in panel (a), the limit cycle dynamics did not evolve along as in the previous examples presented in Sects. 2 to 4. In contrast, the dynamics evolved close to the axis defined by . Panel (a) demonstrates that is a good approximative quantity to characterize the orientation of the limit cycle dynamics in the x–u space. Panel (b) shows , as well. The limit cycle dynamics evolved close to a plane defined by the vector axis , demonstrating again that is a useful characteristic measure.

Fig. 12

State dynamics of the coupled oscillator model defined by Eq. (79). Panels a and b show x(t) (solid line) and u(t) (dotted line) for the initial period of 80s (panel a) and the ”late” period of 80 to 100s (panel b). The figure demonstrates that the model (79) that involves a slow (p) and a fast (s) variable indeed settles down into a limit cycle dynamics. For parameters and initial conditions, see text

Fig. 15

Phase curves, eigenvectors, and effective vector of the coupled oscillator model (79). Panel b: Simulated phase curve u(x) (solid black line), eigenvector axes and (gray dotted lines), and axis (gray dashed line) are shown in the x–u space. Panel b: The solution of Eq. (79) is plotted in the 3D space given by x, u and as phase curve (solid black line). The axis defined by is shown as well (gray dashed line). The figure demonstrates that the orientation of the limit cycle can approximately be described by means of the effective vector and, in doing so, illustrates graphically Eq. (86)

Multifrequency synchronized oscillations on tori

The approach can also be applied to study limit cycles that involve two commensurable frequencies (i.e., frequencies whose ratio can expressed as a rational number like n : m). Limit cycles of this kind may form closed orbits on tori. Let us consider three coupled second-order systems that involve the 3D state vector and evolve likewhere the right-hand side functions (in analogy to Eq. (17)) in general depend on x, u, w and the corresponding first-order derivatives. Let us assume the origin , , , is a fixed point. The model parameters are assumed to depend on a bifurcation parameter, such that at a certain critical value, the fixed point becomes unstable and a Hopf bifurcation takes place. With respect to the fixed point, the evolution equation (87) can (in analogy to Eq. (25)) be cast into the formwhere A, B are matrices. Our approach applies when A and B are simultaneously diagonalizable matrices. In this case, they exhibit the same three eigenvectors that will be denoted by , , and . Moreover, the matrices A, B can be expressed likewith the help of the diagonal matricesand the diagonalization matrix(see also Eqs. (19) and (20). Since (by assumption) exists, the eigenvectors are linearly independent and the state vector as function of time can be expressed in the eigenvector basis likewhere will be referred to as amplitudes again. Let us express the remainder vector likethen (in analogy to the considerations worked out in Sec. 3) from Eq. (88), it follows thatIn summary, in the case of interest the original model (87) can be transformed into Eq. (94). The approach is most powerful if at the Hopf bifurcation point the dynamics of s is characterized by two eigenvalues with negative real parts and is such that s converges to zero. In contrast, the and dynamics is assumed to exhibit eigenvalues with positive real parts.

In the context of the modeling of a n : m multifrequency limit cycle, it then follows that the coupled system , entirely describes the dynamics on the limit cycle. In this scenario, the system at hand exhibits one stable eigenvector and two unstable eigenvectors and and the dynamics on the attractor is given byThe reduced amplitude space is two dimensional.

In the application that will be shown next, for sake of simplicity, it will be assumed that , , and constitute an orthogonal basis. For systems for which this is not the case, a bi-orthogonal basis as described in Sect. 3 can be used.

A class of relative simple multifrequency models, so-called canonical-dissipative (CD) coupled multifrequency oscillators, has been recently introduced [55, 56]. CD models exhibit the advantage that they can be treated to a large extent analytically and include a variety of special cases that have been previously discussed in the literature (for details, see Refs. [55, 56] again). To illustrate by a numerical example the general ideas leading from Eq. (88) to Eq. (95), we used for and the Holt model of the CD framework, which is a 1:2 multifrequency limit cycle oscillator. In addition, we put . The somewhat lengthy equations can be found in the appendix 7.3. Here, we restrict ourselves to discuss the simulation results that are presented in Figs. 16 and 17. We solved the coupled oscillator model defined by Eqs. (88) and (93) (using an Euler forward scheme with time step 0.001s) for the functions and and parameters listed in the appendix. For those parameters, we had , , . The variables x, u, w were assumed to be dimensionless quantities.

Fig. 16

State space (panel a) and amplitude space (panel b) solutions of the Holt 1:2 multifrequency oscillator model based on Eq. (88). Panel a: state dynamics in terms of x, u, w (from top to bottom). Panel b: amplitude dynamics in terms of s(t) (top), (solid line bottom), and (dotted line bottom). The figure illustrates that the dominant dynamics takes place in the 2D space spanned by the unstable amplitudes and . See text for details

Fig. 17

Phase curves, eigenvectors, and limit cycle attractor of the Holt 1:2 multifrequency model. Panel a: The simulated phase curve u(x) (solid black line) and the eigenvector axes and (gray dotted lines) are shown in the x–u space. Panel b: The solution is plotted in the 3D x–u–w space as phase curve (solid black line). Only the final 5s of the simulation are shown. The axes of and (gray dotted lines) and (black dashed line) are depicted as well. The figure demonstrates that the 1:2 limit cycle is located in the 2D space spanned by the axes of the unstable eigenvectors and and, in doing so, illustrates Eq. (95)

Panel (a) of Fig. 16 shows the state variables x, u, w as functions of time. We found that after a transient period, the variables settled down in an oscillatory rhythm that involved two frequencies with a 1 : 2 ratio (as expected). Panel (b) of Fig. 16 shows the amplitudes s, , and as functions of time. As expected, s decayed to zero. In contrast, and converged to a stable multifrequency limit cycle, where oscillated with twice of the frequency of . Figure 17 shows the dynamics in the 2D x–u state space (panel (a)) and the 3D x–u–w state space (panel (b)) and illustrates the role of the eigenvectors. Panel (a) show the phase curve u(x) and the direction specified by . The direction specified by is also shown in the x–u space and (for the set of eigenvectors used in our simulation) happened to be identical to the direction specified by in the x–u subspace ( and only differed with respect to their w-components). As expected, it was found that the dynamics converged towards the direction specified by the unstable vectors and . Panel (b) shows the phase curve of in the 3D space together with the directions specified by all three eigenvectors. Note that the direction related to the stable eigenvector (dashed black line) is perpendicular to the 2D plane specified by the directions of the unstable eigenvectors and (dotted gray lines). In panel (b), only the final 5 s of the dynamics are shown that are assumed to take place approximately on the limit cycle. It was found that the limit cycle (that forms some kind of U or parabola) lies in that 2D plane spanned by the unstable directions. Finally, note that in the 6D space spanned by x, u, w and the corresponding time derivative variables the limit cycle would form a closed orbit that lies on a torus.

Discussion

Unstable eigenvectors and their amplitudes have been proven to be key analysis tools for a variety of studies on nonlinear systems ranging from studies in physics, chemistry, and biology [1, 2, 4, 10, 11, 54, 57–59] to studies of human reactions [13] and research on COVID-19 [18-20]. In the current study, a new application of the concept of unstable eigenvectors and amplitudes was presented for coupled oscillators with simultaneously diagonalizable matrices. In this context, a general class of coupled oscillators was presented with stable limit cycles that are aligned along unstable eigenvectors and whose dynamics on those limit cycles is determined by reduced (one-variable) amplitude descriptions. Three worked-out applications in the field of electric circuits, neuroscience, and cell biology were presented as well.

Coupled oscillatory units such as the Elwakil–Salama circuit (see Sect. 2) have frequently been studied with the help of numerical simulations [28, 29, 60] or by deriving equations for the phase dynamics and discussing phase stability (e.g., Ref. [61, 62]); see also the introduction. In contrast, the proposed eigenvector analysis presented in the current study offers an elegant analytical approach on a level that is ”closer” to the original state space description as compared to the phase dynamics description. In particular, a key benefit of the approach is the one-to-one mapping between the original state space and the amplitude space; see Eq. (13). Features of the original model may be studied in two amplitude sub-spaces that are completely decoupled as in the case of the Elwakil–Salama circuit (Sect. 2) or, alternatively, may be studied with the help of an approximative single-variable description (Sect. 5.1). In general, the amplitude space description may allow for analysis steps that are not readily available on the level of the original state variables. Therefore, the amplitude space description increases the repertoire of tools to analyze electric circuits.

The dynamics of coupled Fitzhugh–Nagumo neurons has frequently been studied by means of numerical methods [31-35] and/or by studying the eigenvalue spectrum [30] and the phase dynamics level [34]. As mentioned in the introduction, these approaches have produced useful insights into possible impacts of neuronal couplings. In particular, in such studies, more comprehensive models have been investigated that exhibit additive terms in Eqs. (40) and (42) and/or in Eqs. (41) and (43). For example, Shim and Husbands [32] considered the asymmetric case for which such additive parameters differ across the two coupled neurons. While it is beyond the scope of the present study to address such more sophisticated features of the coupled Fitzhugh–Nagumo model, future studies may take advantage of the approach presented in Sect. 4.1. In particular, in view of the worked-out example presented in Sect. 4.1, it might be worthwhile to re-examine the cases of the so-called groups C and F presented in the study by Shim and Husbands [32] that seem to exhibit unstable eigenvector axes similar to the ones discussed in Sect. 4.1.

Section 4.2 offers a novel approach to study two coupled Goodwin oscillators. Previous research on coupled Goodwin oscillators focused in particular on the coupling of gene regulatory pathways across cells that are putatively involved in establishing a molecular circadian ”clock” (e.g., Refs. [36, 49]). However, oscillatory signaling pathways have been identified also in other contexts. For example, there is evidence that embryonic development involves limit cycle oscillators on the level of individual cells that cause the period growth of somite segments of zebrafish and the period outgrowth of limbs in chicken [63, 64]. Therefore, the analysis presented in Sec. 4.2 may be applied specifically to study inter-cellular coupled pathways during embryonic development.

The applications presented in Sects. 2 and 4 have been limited to discuss oscillatory systems described by second-order and third-order differential equations. However, in general, and, in particular, in the context of the Goodwin oscillator, fourth-order differential equations [36, 49] and higher order differential equations [65] describing oscillatory units have been used in the literature. Therefore, the current study may open a new avenue to identify the underlying, essential low-dimensional dynamics of such high-dimensional oscillatory models.

Our study focused on limit cycle attractors. In Sect. 5.2, multifrequency m : n limit cycle attractors have been addressed that involve two commensurable frequencies. However, our approach may also be applied to address quasiperiodic oscillations on attractor tori or even chaos. Here, quasiperiodic oscillations are understood as oscillations that involve two incommensurable frequencies (i.e., frequencies and whose ratio corresponds to an irrational number). The phase curve of such an oscillation does not form a closed orbit but typically occupies the entire surface of an appropriately defined torus. While in Sect. 5.2 and the appendix 7.3, the basic CD framework developed in Ref. [55] for the example of a 1:2 polyrhythm has been presented, a more sophisticated and more flexible CD framework has been presented in Ref. [56]. A close inspection of Ref. [56] reveals that this latter framework may be used to demonstrate (just as in Sect. 5.2) the possibility to address quasiperiodic attractors (in the sense defined above), which is, however, beyond the scope of the current work. Similarly, the dynamics of chaotic systems that involve simultaneously diagonalizable matrices may be analyzed using the approach presented in Sects. 2 to 5. To be clear, as discussed at the end of Sect. 3, the approach suggested in the current study does not apply to all kind of chaotic systems. The requirement that the systems of interest exhibit simultaneously diagonalizable matrices imposes certain constraints. However, if (limit cycle, quasiperiodic, or chaotic) systems satisfy those constraints, then the suggested approach can reduce the dimensionality of the problem at hand and help to identify the subspace in which the relevant dynamics takes place. In doing so, new insights might be obtained. In addition, such systems may be used as testbed systems to clarify theoretical concepts or check the performance of numerical methods.

The current theoretical study may also motivate experimental research to test the key predictions that should hold not only for the models that have been discussed in Sects. 2 to 5 but also for similar models that have not been explicitly addressed. In particular, coupled electric units as in Ref. [29] may be examined with respect to their transient dynamics to study the approach towards the relevant unstable direction and with respect to their limit cycle dynamics to fit the dynamics to a single-variable model of an appropriately defined unstable amplitude.