Entrainment of Weakly Coupled Canonical Oscillators with Applications in Gradient Frequency Neural Networks Using Approximating Analytical Methods

Solving phase equations for systems with high degrees of nonlinearities is cumbersome. However, in the case of two coupled canonical oscillators, that is, a reduced model of translated Wilson–Cowan neuronal dynamics, under slowly varying amplitude and rotating wave approximations, we suggested a convenient way to find their average relative phase evolution. This approach enabled us to find an explicit solution for the average relative phase of the two coupled canonical oscillators based on the original neuronal model parameters, and importantly, to find their phase-locking constraint. This methodology is straightforward to implement in any Wilson–Cowan-type coupled oscillators with applications in gradient frequency neural networks (GFNNs).


Introduction
Complex systems transit between various dynamical states. The overall dynamical description of a complex system in many cases can be reduced to a low-dimensional set of nonlinear ordinary differential equations [1]. In cases of a biological complex system, for instance, a nervous system, the descriptive dynamical models could be a Hodgkin-Huxley model that gives a detailed biophysical description of action potential initiation and propagation in neurons, and a Wilson-Cowan neural mass model that describes the dynamics of excitatory and inhibitory neural populations. However, it is rare or even impossible to know precisely all the parameters that are employed in these types of models. To address this issue, one can utilize model reduction techniques to obtain a canonical description of the main models. A canonical model is the simplest class of equivalent dynamical models in analytical terms, and can be derived by methods such as normal form and center manifold theory. For example, the simple spiking model of Izhikevich is a canonical model of Hodgkin-Huxley that could be obtained by either nullclines or current-voltage relations. The canonical model approach is useful for studying the behavior of the entire family of models if the information about most of their parameters is not complete or well known [2].
The canonical model for gradient frequency neural networks (GFNNs) or in short, the canonical model, is a simplified mathematical model that describes the dynamical properties of oscillatory neural populations near a Poincaré-Andronov-Hopf bifurcation [3,4]. Models of cochlear outer hair cells assume critical oscillations based on Hopf bifurcation. For example, it is proposed that in the mammalian cochlea, the basilar membrane is driven by a set of critical oscillators (most likely outer hair cells) and the response of a critical Hopf oscillator has considerable bearing on human psychoacoustics [5]. The canonical model is capable of processing time varying external stimuli, and facilitates the study of nonlinear time-frequency transformations [3], which is important in auditory signal processing. Bifurcation analysis and dynamical properties of a periodically-forced canonical oscillator have been investigated in detail [4]. Such analysis is of importance because the canonical model describes the state transitions of the system between quiescence and spontaneous activity. However, finding the analytical solutions of a canonical model that contains a high amount of nonlinearity is cumbersome. Investigators have introduced and applied approximating analytical methods such as the homotopy perturbation method [6,7], the variational iteration method and more [8]. In this manuscript, we study the average relative phase dynamics of two coupled canonical oscillators while utilizing slowly varying amplitude and rotating wave approximations [9].
Normal form theory, along with center manifold theory, are two approaches that mathematicians employ to simplify dynamical systems. Both of these approaches reduce dimensionality and eliminate the nonlinearities up to the desired order. The method of normal forms can be traced back to the Ph.D. thesis of Henri Poincaré (1929) [10]. It is a powerful tool with which to study the stability and bifurcation analysis of nonlinear differential equations, which employs successive, near-identity nonlinear transformations to obtain a simple form and is restricted to systems which do not have perturbation parameters (unfoldings) [11]. The structure of normal forms is determined entirely by the linear part of the vector field under study. Here we employ conventional normal form theory (CNF), which uses the near-identity change of variables to find the canonical models of single and coupled Wilson-Cowan oscillators.
Our analysis in this paper is categorized under four sections. In Section 2.1 that forms an overview, we start with a Wilson-Cowan oscillator and obtain its normal form. Then we restate the theorems that facilitate deriving a canonical model for single Wilson-Cowan oscillators (the canonical model). In Section 2.2 we extend the methodology of the near-identity change of variables in normal form theory for single oscillators of Section 2.1 to obtain a canonical model for two coupled Wilson-Cowan oscillators. By utilizing polar transformation and decomposing the real and imaginary parts of the system, we obtain radial and angular equations of motion for each of the canonical oscillators. These equations are highly nonlinear and not straightforward to solve analytically. We propose an ansatz in Section 3 that is valid under slowly varying amplitude and rotating wave approximations and leads to a dynamical description of the average relative phase between the two coupled canonical oscillators. This proposed method facilitates working with the cumbersome radial and angular equations derived in Section 2.2. Finally, we solve the average relative phase equation explicitly. This solution is validated by taking its time limit to check if we have correctly obtained the steady state.

Wilson-Cowan Model
The Wilson-Cowan model is a model of neural oscillation based on excitatory and inhibitory neural populations. Neurons are called excitatory if they increase the probability of a postsynaptic action potential occurring, and inhibitory if they decrease this likelihood. According to Hugh D. Wilson and Jack D. Cowan, all nervous processes of any complexity are dependent upon the interaction of excitatory and inhibitory cells [12]. Consider an interconnected pair of excitatory and inhibitory neurons. This pair is called a Wilson-Cowan neural oscillator, and is considered one of the basic mechanisms of the generation of oscillatory activity in the brain [13]. Each unit in the oscillator represents a local population of neurons, in the cerebellum, olfactory cortex and neocortex. By letting x i and y i denote activity of the ith excitatory and inhibitory neurons, respectively, the Wilson-Cowan neural model has the form: where µ x and µ y > 0 are the membrane time constants and τ x and τ y are refractory periods of excitatory and inhibitory neurons, respectively. a ij , b ij , c ij and d ij > 0 for i = j are synaptic coefficients. b ii and c ii are synaptic too, since they denote interactions between excitatory and inhibitory neurons within the ith neural oscillator. a ii and d ii are feedback parameters, which can be positive or negative. Parameters ρ x i and ρ y i denote the external inputs from sensory organs and other regions of the brain to the ith excitatory and inhibitory neurons respectively. S is a sigmoid function that is continuous and monotonously increases. For simplicity and without loss of generality, we consider the Wilson-Cowan neural model in the special case τ x = τ y = 0 and µ x = µ y = 1. That is: We first derive the canonical model for single translated Wilson-Cowan oscillator without external inputs, ρ x = ρ y = 0 in Equation (2):ẋ = −x + tanh(ax − by) Third-order Taylor expansion of the hyperbolic tangent function in Equation (3) leads to the translated Wilson-Cowan model expanded to cubic order: Our strategy in this study (via CNF) is based on the fact that changing the coordinates must not affect the stability of the origin. In fact, we only use coordinate transformations that map the origin to itself.

The Normal Form
Locality, the inclusion of nonlinear transformations and the dominance of the linear part of the system are three requirements of normal form theory [10]. In other words, nonlinear coordinate transformations (near-identity change of variables) are to be performed in the vicinity of a fixed point (in case of Equation (4), the origin), and the structure of this normal form is determined entirely by the linear part of the system under study.
To find the normal form of Equation (4), we map the Cartesian coordinates (x, y) on the complex domain by introducing s = x + iy ands = x − iy. By differentiating s with respect to time and substituting forẋ andẏ given in Equation (4) we get: We have used the fact that x = s+s 2 and y = s−s 2i . Now, we perform a near-identity change of variables, i.e., substitute for s by z according to the transformation: where p and q are positive integers whose sum gives us the highest order of nonlinearity present in the desired canonical model and κ is a complex number. This is an essential step in the normal form theory which gives us the chance to eliminate the desired terms by setting κ appropriately, and consequently simplify the model in terms of nonlinear terms present in the expansion. Finding such values for κ is presented in Appendix A. Rearranging Equation (6) in terms of z yields: If we differentiate this equation with respect to time, we get an expression forż that without considering quintic and higher terms (p + q = 3) is: There are two choices for p and q, either (p, q) = (1, 2) or (p, q) = (2, 1). Both cases are presented in Appendix A, along with detailed calculations. In either case, the general normal form of Wilson-Cowan model (Equation (3)) can be expressed as: with the coefficients presented in Appendix A. We replacedż byż s to indicate that this normal form is for a single oscillator. This CNF can describe all kinds of bifurcations.

The Normal Form about a Hopf Bifurcation Point: The Canonical Model
Here we utilize two theorems from local dynamical systems theory [10] to reduce the normal form in Equation (9) into the canonical model, the model that describes the system about a Hopf bifurcation point.
Note that in the near-identity change of variables we made a transformation of the form z −→ z + h(z,z) where h is third-order in z andz.

Theorem 1.
In transformations of the type z −→ z + h(z,z) for some h = z pzq at the bifurcation point of Hopf type we must not have p − q = −1 [10].
This can never happen if p and q are both even numbers. Therefore, all even-order terms can be removed from the expansion. This theorem guarantees that even if the expansion of the sigmoid function S contains even-order terms, they can be removed for the special case of Hopf bifurcation. This condition also suggests that it works only with p = 2 and q = 1, so Equation (9) can be written as: This is still not the simplest normal form of Hopf bifurcation.
Theorem 2. The cubic terms except z 2z can be eliminated for the Hopf normal form [10].
This is true if and only if: which means a = −d and b = c = ω where ω is a real number. Substituting these results in the other coefficient A 1 leads us to: By introducing α = a − 1 and β = A 21 , Equation (10) becomes: This is the canonical model for a single Wilson-Cowan oscillator without external input in agreement with the results in [13]. In other words, the canonical model (Equation (13)) is the description of a Wilson-Cowan oscillator (Equation (3)) about a Hopf bifurcation point. We compare the nullclines of the translated Wilson-Cowan oscillator with the ones obtained by the canonical model (setż = 0 in Equation (13)) in Figure 1. The nullclines of the models coincide in the vicinity of the origin (the fixed point), which is the point of interest. Here

The Canonical Model of Two Coupled Oscillators
Now we consider two isolated Wilson-Cowan oscillators that are coupled. We write down synaptic and feedback parameters as follows: Additionally, we assume there are no external stimuli to the system, so ρ x 1 = ρ x 2 = 0. Assuming above conditions in Equation (2) leads to: For small synaptic coefficients w 1 , .., w 4 , i.e., to the order of 1, we can rescale the system as w k = c k for k = 1, 2, 3, 4. By choosing tanh as the sigmoidal function S and rescaling, the translated Wilson-Cowan model for two coupled identical oscillators is described by [13]:

The General Normal Form
As we did in Section 2.1, we expand the tanh functions to cubic terms for both oscillators and write the dynamics of each oscillator (i = 1, 2) in the complex domain: The expansion of Equation (16) leads to some terms that describe inherent dynamics of each oscillator (ṡ is ) and other coupling terms (ṡ ic ), which are addressed in details in Appendix B. So in brief: By doing a near-identity change of variables: where m and n are positive integers that satisfy m + n = 3, we obtain the final normal form of these two coupled oscillators:ż Some lengthy but straightforward algebra leads to four cases forż ic depending on the values of p, q, m and n, which is completely presented in Appendix B. The two cases forż is have been already presented in Section 2.1.

Hopf Normal Form: The Canonical Model
As we argued in Section 2.1, Theorem 1, that at the bifurcation point we must not have m − n = −1. Thus, among all four choices of p, q, m and n (presented as four cases in Appendix B), only the case for p = m = 1 and q = n = 2 is acceptable for representing Hopf bifurcation (case (I) in Appendix B). Therefore, the full canonical model that exhibits Hopf bifurcation for two coupled identical Wilson-Cowan oscillators about the fixed point up to quintic terms is: where i = j and all coefficients A 21 , B and A k , B k for k = 1, 2, ..., 8 depend on synaptic and feedback parameters described in the original Wilson-Cowan model in Section 2.1.1 and rephrased in Section 2.2 for the case of two coupled oscillators. All the coefficients have been calculated and are presented in Appendices A and B.

Radial and Angular Equations
Starting with the canonical model (Equation (20)), we can decompose the system to find the dynamical descriptions of each oscillator's amplitude and phase: Considering the polar transformation of z i : z i = re iφ i , represents the canonical model in polar coordinates. We provide all the algebraic calculations in Appendix C and the final coupled nonlinear differential equations that govern the dynamics of the amplitudes and phases are: where i = j and δ is the imaginary part of constant A 21 . We can derive relationships for the relative amplitude and phase of the two oscillators by defining φ(t) = φ 1 (t) − φ 2 (t) and r(t) = r 1 (t) − r 2 (t) as relative phase and amplitude, respectively. We start withṙ =ṙ 1 −ṙ 2 and then substitute for r 1 and r 2 obtained in Equation (22). We do the same for the relative phase equation:φ =φ 1 −φ 2 using Equation (23).
In the special case of r 1 = r 2 the relative amplitude dynamics is described by: where the value of n 5 is given in Appendix C. Solutions of Equation (24) are shown in Figure 2 for different values of the bifurcation parameter α. By setting Equation (24) equal to zero, we found the fixed points, which are either r * = 0 or φ * = 0, π. We demonstrate the solutions of φ for different values of the bifurcation parameter α in Figure 3. For the two cases of α < 0 (damped) and α = 0 (critical) the two oscillators become in-phase with the relative phase of φ = 0. Obviously, the rate of change of φ for the damped oscillation is greater than the critical one. When α > 0 (spontaneous), the two oscillators become anti-phase. In this figure a = α + 1, b = c = ω = 10, d = a − 2, = 0.1, c 1 = c 2 = c 3 = c 4 = , β = A 21 and δ = 0.1. The initial conditions were r 1 = r 2 = 1, φ 1 = π and φ 2 = π/6.

Results
This section contains the main results of this paper. Here, we suggest an analytical solution of Equation (20) using approximating methods. In practice we were interested to find the relative phase dynamics of the two coupled canonical oscillators described above. We considered the compact format of Equation (20) which describes the time evolution of two coupled canonical oscillators (k = 1, 2). Finding the analytical solutions for z is not straight forward. To address this issue, we utilize the ansatz: as an approximate solution of this behavior with c k (t) = r k e iφ k (t) as the time dependent amplitude of oscillations. φ k and ω k are the phase and frequency of the oscillations of each oscillator, respectively. This approximate solution is true under the assumptions of slowly varying amplitude and rotating wave approximations; i.e., the time scale of their dynamics is much larger than 1/ω k (ċ k ω k ), and frequencies bigger than ω k are neglected [1,9,14]. Differentiating this solution with respect to time yields:ż Knowing that:ċ Equation (27) becomes:ż Without loss of generality, we assume equal amplitude for both oscillators, which means r 1 = r 2 = r. Note that even if the amplitudes are not equal, the relative phase dynamics do not get changed since we are analyzing the system under the assumption of slowly varying amplitude approximation; i.e.,ṙ(t) ω. We note here the equations for the first canonical oscillator; the equations for the second oscillator are the same by switching the subindices. The left-hand side of Equation (25) after multiplying by e −iφ 1 e −iω 1 (t) and regrouping is: Equation (30) contains both oscillatory terms (exponential functions) with integer multiples of ω 1 as the frequency and non-oscillating terms that do not contain these exponentials. By integrating over a time period given by T = 2π ω 1 (averaging), all oscillatory terms vanish and only the non-oscillating terms remain: We consider the right-hand side of Equation (25) and set = 1 for now to ease the calculations.
Recollecting the coupling terms on the R.H.S. yields: × [e i(φ 2 +ω 1 (t)) + e −i(φ 2 +ω 1 (t)) ] Introducing new coefficients as: and substituting into the R.H.S. gives: Multiplying the R.H.S. by e −i(φ 1 +ω 1 t) gives: By averaging, only the first term survives and if we consider not necessarily equal to 1: Now by equating the L.H.S. with the R.H.S. and including the second oscillator by letting 1 → 2 and consequently φ → −φ, we obtain: Subtracting second equation from the first one gives: By defining the relative frequency as ω = ω 1 − ω 2 and recalling that sin φ = e iφ −e −iφ 2i : Hence, the average relative phase dynamics is given by: where C = 2M is the total coupling coefficient and a complex number, due to M. We can reinterpret M and consequently C in terms of the initial Wilson-Cowan model parameters (refer to Appendix B): so: Note that the magnitude of the total coupling coefficient is: It can be described as the coupling strength of the two canonical oscillators. This enables us to access a coupling coefficient based on the coefficients in the original Wilson-Cowan model. As one can see, the value of C is independent of c 2 and c 4 . This means that the average phase dynamics of the system do not depend on c 2 and c 4 coefficients. C being a complex number indicates the rotation of the two canonical oscillators in phase space or a delay between their relative oscillations. Integrating Equation (39) gives us more insights: There is a pole at φ = sin −1 ( ω C ). This is the fixed point of the system that gives the steady-state solution of the relative phase. We will confirm it later by taking the limit of φ as time tends to infinity. In other words: Equation (44) gives us the constraint of having convergent solutions; i.e.: This is the convergence interval in which we evaluate the integral of Equation (43). Before solving this integral we consider the special case when the two oscillators have the same frequency, i.e., ω = 0, such that the integration process becomes simpler: where c is the constant of integration and depends on the initial conditions. We solve this equation for φ: hence, the relative phase is: This result can be written in terms of the integration time interval τ = c C : Some plots of the relative phase are shown in Figure 4 for different initial conditions. The relative phase (φ) as a function of time for different initial conditions when C = 0.02. This is a special case when the relative frequency of the oscillators ω is zero, which means they oscillate with the same frequency. One can see that they become phase-locked as time increases. In this special case, since ω = 0, the fixed relative phase will be π.
In the general case of Equation (43) in which the canonical oscillators do not necessarily have the same frequency, i.e., ω = 0, the solution within the convergence region (| ω C | ≤ 1) is given by either: or: depending on the initial conditions φ 1 (0) and φ 2 (0).
Hence, the average relative phase is given by: or: In either case, taking the limit of this average relative phase as time tends to infinity gives us: This angle is the fixed-point or the steady-state solution that we have already mentioned (Equation (44)). As a proof, recall that sin 2x = 2 sin x cos x = 2 tan x cos 2 x = 2 tan x 1 + tan 2 x and if we set x = φ * 2 , then: In Figure 5 we illustrate some solutions of the relative phase for some different initial conditions. Figure 5. The evolution of relative phase when ω 1 = ω 2 . The relative phase (φ) is a function of time for different initial conditions while C = 2 and ω = 1, so we can see the convergence of the solutions and phase-locking phenomenon. In this special case, since ω = 1 and C = 2, the fixed relative phase will be π 6 .
The convergence of the average relative phase of the two canonical oscillators to a fixed point is an indicator of phase-locking. After some amount of time that depends on their initial conditions, the two oscillators become phase-locked; i.e., they tend to oscillate with a fixed relative phase.

Discussion
The canonical model that we studied in this paper is the normal form of the Wilson-Cowan oscillators about a Hopf bifurcation point. The Hopf bifurcation in dynamical systems particularly refers to the development of periodic orbits from an equilibrium or a stable fixed point as an order parameter crosses a critical value [15]-in this paper the bifurcation parameter α (refer to Equations (13), (22) and (23)). This is a point of interest in auditory neuroscience, for example, in modeling approaches that have the potential to link neural dynamics to sound perception [16,17]. Other physical examples include the identification of stable oscillations in nonlinear waves [15]. Nevertheless, if another type of bifurcation is of interest, the normal form could be modified to resemble that particular type of bifurcation. This is because of the fact that normal forms are generic and enable us to remove the nonlinear terms that are not necessary going to be present for our specific purpose.
This work provided both an overview and detailed analysis of the coupled canonical oscillators of Wilson-Cowan-type neural field models, along with a novel methodology to check the interesting solutions. As a recapitulation of previous studies in the literature, we re-derived the canonical models of single and two coupled Wilson-Cowan oscillators under no external stimulus, Equations (13) and (20), respectively. As described in the Methods section, we utilized two theorems in nonlinear dynamical systems that give us constraints on the order of transformations during Hopf bifurcation to find the most compact form of the canonical model. The evolutions of amplitude and phase of each oscillator were found by representing the canonical model in the polar coordinates, Equations (22)-(24), and the solutions were plotted in time for different values of the bifurcation parameter, α, along with the calculation of the fixed points. The damped, critical and spontaneous behaviors were seen as the values of α changes from −1 to 0 and then to 1, respectively (see Figures 2 and 3). These were compatible with the behavior of the coupled oscillators at the Hopf bifurcation, as expected. By integrating knowledge and understanding from the basic mathematics and nonlinear dynamics along with our biophysical understanding of the brain, we presented solutions for the entrainment of two coupled canonical oscillators, which are analyzed in the Results section. These solutions are in line with our previous understanding of coupled nonlinear oscillators; however, a more straightforward method is presented using slowly varying amplitude and rotating wave approximations that helps us obtain approximate analytical solutions of the average relative phases of the two coupled oscillators. The time evolution of the average relative phase is given by Equation (49) for an special case wherein the frequencies of the two oscillators are equal. However, we were able to analytically calculate the average relative phase dynamics in a more general regime wherein the frequencies of the two oscillators were not necessarily equal; see Equations (52) and (53). The results were plotted in Figures 4 and 5 for some different values of the initial conditions. The steady-state of the system, the fixed relative phases of the oscillators are given by Equation (55). The main novel aspect of the proposed method is that it enables investigators to directly relate the parameters in the canonical model to the original oscillatory population model of Wilson-Cowan; for example, see Equation (42). On the other hand, using the ansatz in Equation (26) we overcame the cumbersome task of solving the dynamics of the coupled canonical oscillators that traditionally have been addressed by the method of averaging [15].
Bifurcation analysis and dynamical properties of the canonical model, including phase planes, were investigated in detail in a more general regime of forced canonical oscillators [4] and due to their extensivity, were not presented in this study. The effects of the nonlinear terms in the canonical model, specifically on the resonance patterns, were analyzed and compared with the Wilson-Cowan model in an analytical and simulation study [3], and the analysis of higher order nonlinear terms was out of the scope of this paper.
We considered some assumptions in our analysis to ease the calculations. In all of the study, we considered no external stimulus present in the system. Traditionally, external inputs could be considered within the sigmoidal functions' arguments of the neural field model and consequently would be present in the Taylor expansions of those sigmoid functions. However, the external inputs could be simply considered as direct additive terms as in case of the canonical oscillator, which was investigated in detail in a previous study [4]. For the sake of simplicity, we considered the coupled canonical oscillators to be identical. Nevertheless, one can consider them non-identicall and obtain the same oscillatory terms as in the canonical model but with different coefficients. The average relative phase calculations are done under the assumption of small perturbations; i.e., the time evolution of the canonical oscillators' amplitudes are much smaller than their relative frequency. If we are in a regime in which this assumption is not true anymore, the average relative phases of the two canonical oscillators still could be found numerically but not explicitly.
The results of this study contribute to applications such as neural network modelings that target desired or unwanted resonance in the oscillatory system under investigation. The methodology of deriving the single and coupled canonical models presented here could be generalized for any Wilson-Cowan-type activity rate model and any number of oscillators.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: GFNNs Gradient Frequency Neural Networks CNF Conventional Normal Form

Appendix A
Here are all the calculations needed to continue after Equation (8) and obtain the conventional normal form, Equation (9). In Equation (8) we need to substitute forṡ from Equation (5), but before doing that we rewrite Equation (5) in a simpler way: where the coefficients are: With this substitution forṡ, Equation (8) becomes: Now we transform all s ands on the right-hand side of Equation (A2) according to the near-identity change of variable (Equation (6)) to find the final normal form: By choosing p and q such that their addition is the order of the monomials that we want to remove, we can eliminate many of the monomials and get simpler formats. Since here we are interested in removing monomials with orders higher than cubic, p + q = 3, and there will be two choices for p and q; either (p, q) = (1, 2) or (p, q) = (2, 1). Note that we have made the transformation as z −→ z + h(z,z) for some h = z pzq , where p + q is the order of the term that we want to simplify and p, q ∈ N. We consider these choices of (p, q) as case(I) and case(I I).
• case(I) : (p, q) = (1, 2) In Equation (A3) if we set (p, q) = (1, 2) we get: using the fact that zz = |z| 2 . Applying normal form method successively to eliminate terms that contain z 3 andz 3 we obtain: If we set κ = A 4 2Ā 1 , we can remove one more term which is |z| 2z . Then the normal form becomes: Defining A 21 = A 3 − 2κĀ 2 + A 2κ yields the normal form when (p, q) = (1, 2): • case(I I) : (p, q) = (2, 1) In a similar fashion, but this time with (p, q) = (2, 1) we can eliminate other terms. Again by starting off Equation (A3) and collecting the terms we get: and applying normal form method successively to eliminate terms that contain z 3 andz 3 we get: , the last term on the right-hand side cancels out. By introducing A 12 = A 4 − 2A 2 κ + A 2κ we end up with the final normal form for this particular case: In summary:

Appendix B
Here are all the calculations needed to continue after Equation (19) and obtain the conventional normal form. For succinctness we only do the calculations for the first oscillator; the same calculations are true for the second oscillator by replacing subindex 1 → 2.

(A12)
The first bracket is what we already got for single oscillator in Section 2, and the second bracket is because of the coupling. s index refers to the terms for a single oscillator, and c index refers to coupling terms. We just need to work onṡ c1 , since we have already done the first part (ṡ s1 ) in Appendix A: s c1 = (c 1 + ic 3 )x 2 − (c 2 + ic 4 )y 2 + x 1 y 1 (2abc 1 + 2icdc 3 )x 2 − x 1 y 1 (2abc 2 + 2icdc 4 )y 2 − y 2 1 (c 1 b 2 + ic 3 d 2 )x 2 + y 2 1 (c 2 b 2 + ic 4 d 2 )y 2 − x 2 1 (c 1 a 2 + ic 3 c 2 )x 2 + x 2 1 (c 2 a 2 + ic 4 c 2 )y 2 . (A13) Since x 1 = s 1 +s 1 2 and y 1 = s 1 −s 1 2 , we can replace for x 1 and y 1 in the above equation so we get the right-hand side in terms of s ands: where: and p, q, m and n are positive integers. It is the time now for the near-identity change of variables: In this equation we substitute forṡ 1 from Equation (A14) completely, but for the second and third terms on the right-hand side we do not need a complete substitution, as the other terms are of orders higher than cubic, which are not the matter of interest. Hence: By going back to the original transformation (Equation (A16)) and eliminating all terms of order more than cubic (p + q = m + n = 3), we can obtain the normal form for the coupled part of the oscillations: z 1c = B 1 z 2 + B 5 z 2 1 z 2 + B 7 z 1 z 2z 1 + B 3 z 2z1 2 + B 2z2 + B 6 z 2 1z 2 + B 8 z 1z1z2 + B 4z1 (A20)