A frequency-amplitude coordinator and its optimal energy consumption for biological oscillators

Biorhythm including neuron firing and protein-mRNA interaction are fundamental activities with diffusive effect. Their well-balanced spatiotemporal dynamics are beneficial for healthy sustainability. Therefore, calibrating both anomalous frequency and amplitude of biorhythm prevents physiological dysfunctions or diseases. However, many works were devoted to modulate frequency exclusively whereas amplitude is usually ignored, although both quantities are equally significant for coordinating biological functions and outputs. Especially, a feasible method coordinating the two quantities concurrently and precisely is still lacking. Here, for the first time, we propose a universal approach to design a frequency-amplitude coordinator rigorously via dynamical systems tools. We consider both spatial and temporal information. With a single well-designed coordinator, they can be calibrated to desired levels simultaneously and precisely. The practical usefulness and efficacy of our method are demonstrated in representative neuronal and gene regulatory models. We further reveal its fundamental mechanism and optimal energy consumption providing inspiration for biorhythm regulation in future.

Independent and hybrid coordinations in representative biological models. a Amplitude coordination. The rhythmic biological processes exist in genome-scale transcriptional and translational processes such as protein-microRNA (miRNA) circuitry. The protein concentration in a specific generegulatory network changes throughout the day. The variations under distinct circumstances are sketched by the heterogeneous wave-like pattern (i)-(iii) on the right. The projections of the waves on the rightmost represent their amplitudes. The amplitude is varied under (independent) amplitude coordinations. Meanwhile, the frequencies of these processes are not changed. b Frequency coordination. Periodic biological event is generated due to neuronal excitability exhibiting diverse frequencies under distinct circumstances (right). Under (independent) frequency coordinations, the rates of occurrence can be modulated in both ways (acceleration or deceleration) while the amplitude (intensity of the event) remains invariant. c Hybrid coordination. A biorhythm disorder may include both frequency and amplitude disruption, which further influences internal signal transduction. With hybrid coordination, the two quantities are modulated simultaneously, and they are interchangeable with one another. Compared with pattern (vii), (viii) has a higher frequency and a lower amplitude. An opposite situation is illustrated as pattern (ix). Biological events in the space may occur under different boundary conditions. The amplitude projections shown in the rightmost of (a-c) are typical diagrams for the Robin, Dirichlet, and Neumann boundary condition, respectively. outside the main oscillator are "physically" decoupled by redesigning the oscillator, thus they can be modulated independently in extended ranges.
Actually, the regulations of frequency and amplitude of the components inside an oscillator are also significant. Moreover, besides the network topology, the intensities of endogenous interactions affect frequency and amplitude decisively, but there are few related works focusing on this direction 33 . Frequency and amplitude can be coordinated by making a slight intervention on these interactions. This can be realized either internally by varying the system parameters (e.g., degradation rate, synthesis rate, strength of a stimulus, etc.) or externally by a feedback controller (e.g., computer-based microfluidic devices 34,35 ). More importantly, the underlying mechanisms for coordinations as well as their energy consumptions also depend highly on these interactions, which have not yet been uncovered. In this article, we present a universal computational framework for precisely designing such an intervention coordinator. We will try to answer the following questions: What specific form of an intervention should be made (i.e., enhancing or mitigating the original interactions)? How is its intensity and energy consumption? Additionally, with our well-designed coordinator, both frequency and amplitude can be regulated precisely and concurrently (Fig. 1c).
To get a deeper insight into the biological oscillators, dynamical systems techniques based on computational and mathematical modeling are usually applied 21,[36][37][38][39][40][41] . For systems in small scales, spatial movements including molecular and ion diffusion cannot be neglected when modeling because they would cause heterogeneity patterns (Fig. 1a, b). By dividing space into several discrete parts, we can use agent-or compartment-based models (i.e., ordinary differential equations) to simulate such a system 42 . However, these models still have their limitations. We may oversimplify a system if the spatial division is coarse, and the large number of components would increase computational consumption hugely. Contrarily, partial differential equation (PDE) is a more appropriate model to describe systems with variations in both continuous space and time. Especially, reaction-diffusion (R-D) systems are commonly used to investigate dynamics at the molecular, cellular, and tissue levels [43][44][45][46][47] . It is also well known that Alan Turing used both discrete and continuous models in his seminal work on morphogenesis 48 . Moreover, R-D systems can describe simple but significant random walk processes 49 . The Hopf bifurcation plays a pivotal role in appearance of a rhythmic oscillation. It refers to a transition from a quiescent state to an oscillatory state [50][51][52][53] . In R-D systems, it yields spatial, periodic oscillations 54,55 . Owing to diffusive effects and diverse boundary conditions, such rhythmic oscillations may be spatially non-homogeneous (Fig. 1a, b), for which frequency and amplitude coordinations still remain to be explored. It is then reasonable and significant to start our study with the Hopf bifurcation in R-D systems.
In this work, we computationally design a universal coordinator intervening linear interactions for modulating both frequency and amplitude. It is feasible for oscillations arising from the Hopf bifurcation in generic two-component R-D systems with the conventional boundary conditions [e.g., the Neumann boundary condition (NBC), the Dirichlet boundary condition (DBC), and the Robin boundary condition (RBC), see Supplementary Note 1 for details]. The basic idea is to extract and normalize the frequency and amplitude by utilizing the center manifold and normal form theories 54,56,57 and to link them with the intensities of linear interactions. In such a way, we are able to design a useful coordinator leveraging its linear regime [e.g., the Michaelis-Menten (MM) function]. In two representative biological models (Fig. 1a, b), a gene-regulatory network and a neuronal model, we demonstrate the efficacy and practical usefulness of our framework. In light of the estimations computed via our approach, oscillations far from the quiescent state can be modulated as well. The underlying mechanisms and energy consumptions are also discussed, which might be helpful for biological regulation in future.

Results
Coordinating a "cancer network" by MM regulations. As a first example, we consider the cyclic dynamics in a gene-regulatory network involving a microRNA (miRNA) cluster and a protein module 58 . It is called a "cancer network" because the miRNA behaves as an oncogene or tumor suppressor depending on the protein concentration. Its mathematical model abstracts the interaction among the transcription factors and a miRNA cluster 58 whose dimensionless diffusive model is given as (see Supplementary Note 8 for more details on the model) where ϕ = ϕ(x, t) and μ = μ(x, t) represent the dimensionless level of protein and miRNA, respectively. Here we include two diffusive terms (second-order differentiation with respect to the spatial variable x) into the equations to describe the unbiased molecular diffusion. With appropriate parameters (Supplementary Table 3), the model possesses a constant quiescent state (ϕ 0 , μ 0 ) that undergoes the Hopf bifurcation at ϵ = ϵ* yielding rhythmic oscillation when ϵ < ϵ* (Fig. 2). To achieve frequency and amplitude coordinations, we add two terms F 1 (ϕ, μ) and F 2 (ϕ, μ) into the first and second equation to regulate the dynamics of protein and miRNA, respectively. They act as instantaneous regulations on ϕ and μ as follows F 1 ðϕ; μÞ ¼ f 11 Mðϕ À ϕ 0 ; K 11 Þ þ f 12 Mðμ À μ 0 ; K 12 Þ; where M(x, K) = x/(x + K) is a MM regulation and f ij indicate the intensities. The MM regulation was recently shown to be practical and useful in implementing the proportional control in biological molecules 59 . Here we take into account all possible regulations: two self-regulations (ϕ → ϕ and μ → μ) and two cross-regulations (ϕ → μ and μ → ϕ). We need to point out that the regulations in Eq. (2) incorporate the coordinate translations because the dynamics considered in our case oscillates around the quiescent state (ϕ 0 , μ 0 ) and we desire to leverage the information on their displacement from the quiescent state. Actually, such a translated regulation can be divided into two typical MM functions with K ¼ K À x 0 . For the given Michaelis constants (K ij = 5 in our case), we select appropriate regulation intensities f ij to coordinate independently the frequency or amplitude of the oscillation near the Hopf bifurcation (ϵ = 0.08 < ϵ*). By varying the intensities at distinct time, we eventually accelerate the oscillation with a doubled frequency. Concurrently, both amplitudes of protein and miRNA concentrations are kept near constants (Fig. 2a-d). Moreover, with other well-designed intensities, the independent amplitude coordination is also successfully performed (Supplementary Fig. 3 and Movie 5).
Depending on the levels of protein module ϕ, the miRNA cluster is classified as oncogenes or tumor suppressors 58 . When ϕ approximately lies between 2.75 and 3.75, the region is labeled as a cancer zone where hyper-proliferation occurs. Accordingly, the probability of oncogenesis is increased. We find that the oscillation stays away from the cancer zone (i.e., ϕ < 2.75) if it is close to the quiescent state ( Fig. 2a, b). But, when it keeps growing, its amplitude increases and the peak enters the cancer zone. For instance, the peak of the oscillation (when ϵ = 0.05 without a coordinator) shown in Fig. 2e, f (t < 50) lies in the cancer zone. To prevent the oscillation from the entry of the cancer zone, it is possible to use our coordinator with appropriate intensities f ij to suppress its amplitude. For this purpose, we design appropriate coordinators and apply them at different time to decrease the amplitude gradually (50 < t < 300 in Fig. 2e, f). Eventually, the protein concentration stays away from the cancer zone. Thus, the miRNA cluster is no longer classified as oncogenes. Evidently, our designed coordinator is feasible for coordinating both frequency and amplitude. Let us now introduce its theoretical background and a universal computational framework for determining the intensities. : ð4Þ Such a system possesses two self-interactions and two crossinteractions, all of which can be divided into two parts: the linear interactions Au and the nonlinear ones g. The matrix A is associated with the Jacobian matrix and can be acquired in any computational model 21 . It represents the endogenous linear interactions, which consists of four coefficients a ij (i, j = 1, 2) representing associated interaction (u j to u i ). The sign of a ij indicates the role played by the component j (activator or inhibitor), see below. Their magnitudes represent the intensities that are determined by system parameters, such as degradation rate, synthesis rate, Michaelis constant, etc. These parameters are usually different from system to system (see Supplementary Note 9 for more details on the linear interactions). For a given biological system, the endogenous linear interactions always exist in its computational model and have dominant  impacts on the behavior of the oscillation near the quiescent state. Specifically, the intensities of these interactions determine the frequency and amplitude decisively. When translating the coordinated system, our proposed MM regulations [Eq. (2)] are also shifted to a regular MM function. For instance, M(ϕ − ϕ 0 , K 11 ) is translated to M(Φ, K 11 ) with Φ = ϕ − ϕ 0 representing the displacement of the oscillation from the quiescent state. Near the Hopf bifurcation, the quantity Φ is relatively small. The MM function is therefore in its linear regime, that is,

Endogenous linear interactions and their interventions.
Consequently, near the Hopf bifurcation, the introduced MM regulation is analogous to a linear proportional control making interventions on the original linear interactions, thereby our proposed coordinator is feasible for coordinating the frequency and amplitude.
A universal design policy for the coordinator. To design a feasible coordinator, we need to know how the intensities of the MM regulations relate to the frequency and amplitude. Here we apply the center manifold and the normal form theories to accomplish the task (Fig. 3). In the following computational framework, we consider the linear regime of the MM function where the feedback coordinator for Eq. (4) is written as Note that the intensities f ij in Eqs. (6) and (2) (6) is also suitable for other nonlinear functions possessing a linear regime, such as the sinusoidal signal/current and other generalized MM functions We now present the computational framework that includes the following steps. We first perform a persistence analysis (see "Methods") before moving into detailed computations. It guarantees the existence of the oscillation after intervention. This gives us Theorem 1, which provides the first criteria (f 22 = − f 11 ) for designing the coordinator. It indicates that the two self-interactions of a feasible coordinator must be opposite to one another. Thus, a negative feedback loop is established providing the potential to produce rhythmic oscillation in a biological system 29 , reasoning that it ensures the occurrence of the Hopf bifurcation leading to a cyclic oscillation. In practical applications, a stable oscillation is of great significance and interest. Therefore, we also perform a stability analysis to guarantee that the oscillation is still stable after implementing the coordinator (see "Methods").
The second step is to link the intensities f ij with the frequency and amplitude. Therefore, we exploit the normal form of the Hopf bifurcation 57 utilizing the center manifold and the normal form theories. For a periodic oscillation, the Poincaré normal form is always found as (see "Methods" for more details) where λ = μ + iω is the complex eigenvalue corresponding to the Hopf bifurcation and η is the normal form coefficient whose real part, denoted by χ, is the first Lyapunov coefficient. It is the simplest form preserving the information and is usually applied to study qualitative behaviors. Here we extract the quantitative information (frequency and amplitude) from the normal form. It serves as a bridge from the original oscillation to the coordinated one. Therefore, those information are helpful for designing the coordinator. The quantifying procedure is illustrated in Fig. 3. The third step is to modulate both components (u 1 and u 2 ) consistently. We notice that they always share the same frequency, whereas their amplitudes could be distinct. As a consequence, the variations made by the coordinator may be inconsistent, in the sense that the amplitude of u 1 (e.g., protein concentration) is increased or decreased while that of u 2 (e.g., miRNA concentration) is invariant or vice versa. Such a result may be unexpected for a given biological system. To avoid this Fig. 3 The process for quantifying frequency and amplitude of a periodic oscillation. A typical periodic oscillation with two components (u 1 and u 2 ) in a R-D system with non-uniform distribution is shown in the leftmost panel. The diffusive (heterogeneity) effect along the spatial variable x can be clearly seen. By applying the center manifold theory, such an effect is eliminated (see Supplementary Note 5 for details). Then the restricted periodic oscillation is spatially homogeneous (invariant along x-axis) and finally can be projected into a complex plane which is of dimension two, as shown in the second panel. By applying the normal form theory, the periodic orbit is converted to a circle (the third panel), whose frequency [ω/(2π)] and amplitude ( circumstance, we analyze both components when computing the normal form (see "Methods"). We surprisingly find that the modulations on both amplitudes are commensurate with one another if we have: This becomes the second criteria for the design policy. Reasonably, it incorporates the endogenous cross-interactions and their interventions. The ratio of two amplitudes is affected by the crossintensities. To keep this ratio unvarying, the intensities of the interventions must also follow a fixed ratio determined by the intrinsic ones (i.e., a 21 /a 12 ).
In the final step, we derive two algebraic equations for the coordinator. From Eq. (9), the frequency and amplitude are extracted and approximated as ω/(2π) and 2 ffiffiffiffiffiffiffiffiffiffiffi Àμ=χ p , respectively. Note that μ is independent of the intervention coefficients, whereas ω and χ are expressed in terms of f 11 and f 12 , i.e., ω = ω(f 11 , f 12 ) and χ = χ(f 11 , f 12 ). Before going any further, it is worth pointing out that the approximated amplitude is independent of x (i.e., the spatial variable). In spite of this, the coordinator designed here still works for the entire space, because the spatially heterogeneous effect is eliminated during the computation (see Supplementary Notes 5 and 6 for more details) and it does not affect the final coordination.
We then denote by ω 0 and χ 0 the two values for the original oscillation [i.e., ω(0, 0) and χ(0, 0)], and ω c and χ c the coordinated ones. To accomplish the coordination at the desired frequency and amplitude, we only need to solve two algebraic equations: where r F and r A indicate, respectively, the fold changes of the frequency and amplitude with regard to the original quantities. For instance, to coordinate frequency independently, we set r A = 1 and an unrestricted r F . Analogously, we coordinate amplitude independently by setting r F = 1 and a free r A . The frequency and amplitude are coordinated simultaneously if neither of r F and r A is one. Generically, the equations afford an accurate prediction for designing the coordinator ( Supplementary Fig. 5). For higher accuracy, we numerically search more appropriate values close to the predicted ones.
Following the above procedure, we design feasible coordinators for the independent frequency and amplitude coordinations in the "cancer network" (see Supplementary Note 8 for computational details). As mentioned before, Eq. (6) provides accurate estimation of Eq. (2). To verify this fact, we also compare the efficacy of the linear coordinator and that of the MM regulations. Our proposed nonlinear coordinator is indeed well approximated by the linear coordinator even if the oscillation is relatively far from the quiescent state ( Supplementary Fig. 4).
Independent coordination in a neuronal model. To further demonstrate the efficacy of our framework, we also apply it to the FitzHugh-Nagumo (F-N) system, a representative model simulating biological neuron 60-62 (see "Methods"). The parameters we used under distinct boundary conditions and other information are provided in Supplementary Table 5. Following the procedure introduced before, we design some feasible linear coordinators such that the frequency of a periodic oscillation is gradually increased or decreased under different boundary conditions (Fig. 4a,  Given a fixed fold change r F , some computed intervention intensities should be abandoned as their magnitudes are too large for practical applications. Besides, to guarantee that the coordinated periodic oscillation is still stable, we require the coordinator to satisfy two stability conditions: the negativity of the first Lyapunov coefficient χ and the positivity of another measure index c 0 (k i ) (see Theorem 2 in "Methods"). Accordingly, a fixed r F yields two coordinators (f 11 , f 12 ), where f 11 are identical and the two values of f 12 are symmetric to a constant.
We present in Fig. 4b-d the relations between the intervention intensities and the ratio of frequencies. The two coefficients show an almost linear relation. When ω 0 /ω c < 1 (i.e., r F > 1), the large magnitudes of f 11 and f 12 imply that a more intense coordinator is needed for acquiring a higher frequency. As a comparison, it is much easier to decrease the frequency in the sense that the required interventions are relatively moderate.
Analogously, we also accomplish the independent amplitude coordinations with the proposed approach. Under different policies, the oscillation is gradually amplified (for both components) while the frequency remains unchanged (Fig. 5a, b). See Supplementary Movie 4 and Figs. 13, 14, and 16 for more instances. We find that the amplitude fold change r A has a lower bound (e.g., about 0.723 for the NBC case), below which the oscillation cannot be coordinated successfully (Fig. 5c, d). In other words, the periodic oscillation cannot be suppressed to an arbitrarily small size. This is caused by the instability of the oscillation. Although we find some coordinators for every r A < 1, they may violate Theorem 2 (see "Methods"). Consequently, the oscillation becomes unstable. To analyze the lower bound, we investigate the minimum of the index c 0 (k i ), whose negativity implies instability (see Supplementary Fig. 15). In the unstable region (Fig. 5c, d), the oscillation cannot be observed since it would converge to a different state (possibly a steady state) or undergo a finite time blowup.
As claimed before, one feature of our method is the capability to modulate the amplitude in the entire spatial domain by almost identical fold change r A . Specifically, we illustrate in Fig. 5e, f an adequate example. Apparently, for both components (V and W), the amplitudes are increased for all spatial variable x (the vertical axis) indicating that the spatially non-homogeneous property does not influence the effect of our coordinator.
The optimal energy consumption of a coordinator. In practice, the consumption of a coordinating policy for sustaining an oscillation is of great significance. Therefore, an index E is always introduced to measure the performance of a specific control policy 63 , see "Methods" for the one used in our problem. Figure 4e shows the energy consumption of independent frequency coordination in the F-N system with NBC. It is clear that there is a big difference between the two policies. Usually, the blue one is the optimal one that we may use because its consumption is lower. It is worth pointing out that, throughout this work, we only illustrate the coordinations attained by the optimal coordinator. A clearer energy consumption of the optimal frequency coordinator for the F-N model is provided in Supplementary Fig. 11 from which we deduce that the coordinator consumes more energy for acquiring a decreased frequency (i.e., ω 0 /ω c > 1). In Fig. 5d, the energy consumption for the optimal independent amplitude coordination is also given. Apparently, the magnitudes of the index E are exceedingly different for the two cases in the sense that the amplitude coordination consumes more energy than the frequency coordination. Moreover, increasing and decreasing the amplitude yield analogous energy consumptions. In Fig. 6c, the optimal energy consumptions for the four cases are indicated according to their magnitudes.
Coordinating the "bigger" oscillation. Until now, our framework are verified on distinct oscillations in the two models, most of which are the oscillations near the Hopf bifurcation with relatively small amplitudes. Our method, based on dynamical systems tools, provides a reliable design policy for coordinating such oscillations. In practice, an oscillation far from the quiescent state with greater amplitude is, however, more likely to exist and needs to be coordinated. For such an oscillation, our approach provides an estimation of a feasible coordinator. Mostly, it may not be accurate enough to obtain desired frequency or amplitude because they are also affected by the nonlinear terms as the oscillation grows bigger. But, we can vary the four intensities f ij and search a more appropriate combination numerically in the neighborhood of our estimations. Compared with seeking a feasible coordinator in the entire space R 4 , the numerical investigation based on proposed approach significantly reduces the computational consumption. As two examples, we perform the independent coordinations on "bigger" oscillations in both "cancer network" (Fig. 2e, f) and F-N system (Supplementary Fig. 16).
Independent coordinator is an amplifier or a damper. Though we have successfully performed independent coordinations under distinct circumstances, there is still no overview picture on the coordinator configurations. That is, it is not apparent whether the four interventions in Eq. (6) are positive or negative, whereas this information could be significant for understanding the underlying mechanisms. For this purpose, we ignore the exact values of the intensities and focus only on their signs, positive (activator) or negative (inhibitor) (Fig. 6). We observe that, for both models, in order to acquire an increased frequency or amplitude, the required intervention coefficients follow the endogenous ones in the same direction (i.e., f ij ⋅ a ij > 0), while we have f ij ⋅ a ij < 0 for a decreased frequency or amplitude. Specifically, the designed coordinator is a reversible version of the endogenous interactions for suppressing either frequency or amplitude. In summary, the designed independent coordinator is just like a two-way dial, which can be turned up or down. When serving as an amplifier, it The relation between f 11 and f 12 of an optimal coordinator (the one with minimum energy consumption). The solid curve and circles represent theoretical and numerical results, respectively. The oscillation is accelerated (respectively, decelerated) when f 11 > 0 (respectively, f 11 < 0). When f 11 approaches the critical value (approximately −0.25), the period of the oscillation tends to infinity. The oscillation eventually disappears in the darkest region. The coordinating policies applied in a as time increases are highlighted (from left to right) with light blue circles. c, d The values of f 11 and f 12 for obtaining different frequencies. Two distinct strategies are shown in different colors. Note that the horizontal axis indicates the value of ω 0 /ω c (i.e., 1/r F ) for readability. e The average energy consumed by a coordinator to sustain an oscillation at desired frequency. The two colors correspond to the two policies shown in d. f Independent frequency coordination (threefold decrease) in the F-N neuronal model with DBC. The frequency is gradually and independently decreased by varying the coordinator at each time labels t = 200, 400, 600, 800, and 1000. accelerates and reinforces the reaction in a system raising the frequency or amplitude. Otherwise, the coordinator serves as a damper slowing down or mitigating the reaction (akin to breaks of a car).
Hybrid coordination and energy transfer. Besides modulating frequency or amplitude of a biological oscillator independently, sometimes a hybrid coordination is also needed to meet a certain demand (e.g., signal transduction). In this way, the frequency and amplitude can be interchanged with one another. We can design a coordinator for any values of r F and r A as long as the two algebraic equations have a solution and the coordinated oscillation remains stable. Therefore, the hybrid coordination can be easily accomplished by our approach. Figure 7d shows the region of hybrid coordination bounded by the two curves for the independent coordinations (the neuronal model with zero flux at the boundary). Note that the configurations for independent frequency and amplitude coordinations (Fig. 7f, g) obey the mechanisms presented in Fig. 6. For an oscillation moving in the hybrid region, its frequency and amplitude are coordinated in a reverse way. An example for the hybrid coordination is shown in Fig. 7a-c. The oscillation in the beginning (t < 200) has already experienced a frequency coordination. It has the same amplitude and higher frequency (r F = 5/3) as the original oscillation. As time goes on, it undergoes a hybrid coordination along the green curve (circles from left to right) in Fig. 7e. With such a sequence of coordinations, the amplitude is increased (up to r A = 5/3), whereas the frequency is gradually decreased until it reaches the same amount as the original one (1000 < t < 1200). We also computed the energy consumed by the hybrid coordinators (Fig. 7h, i). Along the hybrid coordination, more energy are required from sustaining higher frequency to higher amplitude, because an oscillation with higher amplitude consumes more energy than the one with higher frequency. If the hybrid coordination is performed in an opposite way, some energy can be released. This finding is in accordance with the results shown in Fig. 6c.

Discussion
Over the past decades, it was widely recognized that frequency and amplitude variations are fundamental processes in biological oscillators. While their precise coordinations via either internal regulation or external intervention have great benefits for health, the underlying mechanisms and the coordinating strategies are rarely studied. Especially, there are few works focusing on concurrent regulation of both quantities rather than a single one and  Supplementary Fig. 12 for the time course of W. b The time course at x* = π/2. c The relation between f 11 and f 12 of an optimal coordinator for the independent amplitude coordination. The circles (from left to right) represent the policies applied in a as time goes on. When f 11 is greater (respectively, less) than zero, the amplitude of the oscillation is increased (respectively, decreased). A critical value exists (approximately r A = 0.723), below which the oscillation turns out to be unstable. d The average energy consumed by a coordinator to sustain an oscillation at desired amplitude. When suppressing an oscillation (the light gray region), the energy consumption grows rapidly as the amplitude becomes smaller. Contrarily, the energy growth is not significant for magnifying the amplitude (white region). the intensities of the motif are usually ignored. Through this work, we accomplish this task with a well-designed coordinator intervening the linear interactions. We find that the interventions of linear interactions play dominant roles for the coordinations on both frequency and amplitude. Also, with two concrete and successful examples, we demonstrate that the designed coordinator are heavily dependent on the endogenous linear interactions. According to required outcomes, it is considered as either an amplifier or a damper (Fig. 6). This unprecedented finding may lead us to a better understanding of the underlying mechanism of biological oscillations. Of course, it is a challenge to address the questions: Is our discovery a universal principle? Is there any other mechanism behind frequency and amplitude coordinations? They are beyond the scope of the present article, which could be important directions for future study. The proposed coordinating strategy can be regarded as either an internal or an external intervention. Therefore, our discovery provides significant information for not only the steering of traditional biological models but also the designing of flexible artificial biosystems in future.
Compared with the recent work focusing on the re-design of specific oscillators 32 , our approach pays attention to the oscillator itself and the intensities of its linear interactions. On the one hand, we focus on coordinating the frequency and amplitude of the components inside the oscillator rather than using a main oscillator to control an outside gene 32 . Therefore, the two works have practical usefulness under distinct circumstances. If the main oscillator is more important and needs to be modulated, then our framework is appropriate to accomplish the task. On the other hand, our designed coordinator is able to precisely acquire desirable frequency and amplitude. It provides the possibility of precise coordinations. Though it has not been verified experimentally in this work, it may be realized in future by a welldesigned computer-based microfluidic device. For various biological systems, such as the neuronal model and the gene-regulatory network considered in this work, the endogenous interactions are always extractable, if its corresponding computational model is given (or discovered via machine leaning from data 64,65 ). Therefore, the proposed strategy can be applied to those systems straightforwardly. This work serves as a necessary foundation for precise biorhythm regulation. It may fill the gap between experimental data and theoretical coordinating strategies making a potential contribution for precision medicine.
In the field of dynamical systems, the center manifold and the normal form theories are usually used to carry out qualitative analysis, specifically, to investigate the direction of the Hopf bifurcation and the stability of associated periodic oscillation 57 .
Here we show that they are also powerful for coordinating biological oscillators. We exploit their roles for analyzing quantitative information of rhythmic oscillations. In such a way, the Fig. 6 The underlying mechanisms and energy consumptions of independent frequency and amplitude coordinations. a The optimal coordinator [the one with minimum energy consumption E, see "Methods" for its definition] for independent frequency and amplitude coordination follows the two rules. The four linear interactions a ij (arrow: activation; I-shaped: inhibition) between two components (blue and red circles) are intervened by the linear regime of the coordinator [Eq. (6)] with intensities f ij . For both independent frequency and amplitude coordination, all the four interactions are either enhanced (thicker curves) or mitigated (thinner curves) simultaneously. Thus, the coordinator can be viewed as either an amplifier or a damper. b The nonintervened oscillation in the neuronal model Eq. (14) (blue: V; red: W). All four configurations for independent frequency and amplitude coordination are tabulated below. Only the signs of the linear interactions and interventions are taken into account. If f ij (or a ij ) is positive, then component j can be regarded as an activator to i, otherwise, an inhibitor. It can be observed that all the intervention intensities follow the same sign as that of the corresponding endogenous interaction, if an increased frequency or amplitude is desired (the case of enhancement in a). On the contrary, for acquiring a decreased frequency or amplitude, the designed intervention exhibits as a reversible version of the endogenous one (the case of mitigation in a). c The four intervened examples for each case given in b. The energy consumption for each case is sketched by arrows on the right of the motif. A higher and darker arrow suggests a greater energy consumption. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-26182-2 ARTICLE intervention intensities and desired frequency or amplitude are connected, whereby simple criteria are derived for designing coordinators. Consequently, we can alter the frequency and amplitude independently or simultaneously. With these tools, we also show that the diffusive effects are preserved from the original oscillation during the coordination so that our approach is applicable for the entire spatial domain. In fact, biological processes following certain rules are dynamical systems described by either discrete or continuous mathematical models. Besides the Hopf bifurcation yielding rhythmic oscillation, there are other bifurcations accounting for various biological phenomena 53 .
Analogously, dynamical systems tools can be applied for certain biological coordinating tasks, such as controlling spatial patterns induced by Turing instability 66 and excitability related to canard phenomenon in biological systems with multiple time scales 67 , to name a few.
The present approach sheds some light upon the tunability of heterogeneity biological oscillators and their energy consumption. Though we only consider two-component systems in this work, our systematic approach can be readily adapted to the one with more components following the idea that we can decouple and extract the information of frequency and amplitude via the center The relations between f 11 and f 12 of an optimal coordinator for the independent frequency and amplitude coordination are shown, respectively, by the blue and yellow curves. The gray area between the two curves is labeled as hybrid coordination region, where the frequency and amplitude can be coordinated and interchangeable with one another. e A zoom-in view at the red rectangle in d. The oscillation shown in a is coordinated along the green curve. The policies applied as time goes on are highlighted by circles from left to right. f, g The optimal configurations for the four cases tabulated in Fig. 6b (the curves for ij = 11 and ij = 22 are not visually distinguishable). For an increased (respectively, decreased) frequency (blue) or amplitude (yellow), we have a ij ⋅ f ij > 0 (respectively, <0). h, i The frequency (blue) and amplitude (yellow) variation (r F and r A ) together with the average energy (green) consumed by hybrid coordination along the green curve in d.
manifold and the normal form theories. For instance, the Kim-Forger model describing the circadian rhythm regulated by BMAL1:CLOCK transcription factor and PER:CRY complex 53 can be investigated with the present approach. In this model, the cyclic oscillation also arises from the Hopf bifurcation. More components imply more intervention coefficients. Therefore, we would have more choices to design coordinators for such systems. More importantly, the energy consumption may be further reduced or follows a different rule compared to the one in Fig. 6. Remember that the amplitude coordination is restricted in some range due to the instability of the oscillation. With diverse coordinating policies, this range may also be extended.
Although we show that the linear regime of a well-designed coordinator works well for both "small" and "big" oscillations in the two models, its efficacy may be less than satisfactory for systems with greater nonlinear coupling effect. Also, a coordinator sometimes cannot be instantly implemented due to delay. Therefore, besides the linear interventions, we can also take advantages of nonlinear interventions (e.g., the Hill function of multimer) or the one with time lag to improve our results in future. On the other hand, an oscillator over a two-or threedimensional spatial domain could be considered and some boundary control techniques 68 are also possible for coordinating such oscillations.

Methods
Generic R-D model. The periodic oscillation arising from the Hopf bifurcation are studied in this work. Thus, for every R-D system investigated here, it possesses a constant stationary solution (i.e., homogeneous quiescent state). By appropriate transformation and Taylor expansion, this solution is moved to the origin and the system is written as Eq. (4).
is the nonzero diagonal matrix consisting of non-negative diffusive coefficients, A(ϵ) is an n × n real matrix, and g(u, ϵ) includes nonlinear terms that satisfies g(0, ϵ) = 0. Particularly, we focus on the case of n, the number of components, as 2 throughout the work. The methods proposed in this work can be extended straightforwardly to a system with more components or spatial variables.
Coordinator configuration. The coordinator used in the computational framework is Eq. (6). For simplicity, we define a linear operator as LðϵÞ :¼ DðϵÞ∂ 2 =∂x 2 þ AðϵÞ þ F. The coordinated system is then written concisely as Eigenvalues of the operator L. For all the three boundary conditions considered in this work, the Laplacian operator in L possesses countably infinitely many eigenvalues k 2 i ðk i ≥ 0; i 2 N 0 Þ, and they are always ordered as 0 ≤ k 2 0 <k 2 1 <k 2 2 < Á Á Á 69 . Further, for each k i , we deduce the eigenvalue problem and ψ k i is the eigenvector corresponding to λ k i . Then, every λ k i is also an eigenvalue of L and can be solved from the characteristic equation given by For the details on the corresponding eigenvector, see Supplementary Note 2.
Persistence and stability analysis. Depending on the particular parameters and boundary conditions, the periodic oscillation arising from the Hopf bifurcation would be either stable or unstable 57 . First, to guarantee its existence with or without a coordinator at the same bifurcation parameter, we have the following proposition (see Supplementary Note 3 for more details).
As the Hopf bifurcation occurs at the same critical value for any F, we deduce the theorem below (see Supplementary Note 4 for a detailed proof). Theorem 1. The Hopf bifurcation of u ≡ 0 in Eq. (11) always occurs at the same system parameters if and only if the matrix F satisfying f 11 + f 22 = 0.
Second, to guarantee that the Hopf bifurcation yields a stable oscillation, we also have the following proposition (see Supplementary Note 3 for more details).
Proposition 2. When ϵ = ϵ*, the operator L has a unique pair of purely imaginary eigenvalues yielding the Hopf bifurcation. All other eigenvalues strictly lie in the left-half complex plane.
We now analyze the spectrum of L and determine ϵ*. According to Proposition 2, the roots of Eq. (12) must be non-zero and their real part are non-positive (for all k i ). Then, a simple manipulation (Supplementary Note 4) implies the following theorem.
Theorem 2. c 0 (k i ) is strictly positive for all k i .
Generically, the quiescent state u ≡ 0 undergoes infinitely many Hopf bifurcations. That is, for each k i , there is a critical value such that Eq. (12) has a pair of purely imaginary roots. Consequently, there exists a sequence of periodic oscillations arising from each Hopf bifurcation. It is then natural to raise the question: Which oscillation should we coordinate its frequency and amplitude? Thanks to the following theorem (proved in Supplementary Note 4), the only case that we need to investigate is i = 0. All other Hopf bifurcations yield unstable oscillations that is of less interest in practice. We finally have the two theorems below.  Normal form and quantitative information. Having analyzed the stability of periodic oscillation, we then compute its normal form, from which its frequency and amplitude are extracted. We do this by applying the projection method 57,70 . According to Theorem 3, the pair of purely imaginary eigenvalues is solved from Eq. (12) when i = 0. For readability, we drop the subscript and represent the two eigenvalues as λ and λ. Then, for ϵ sufficiently close to ϵ*, they can be written as λðϵÞ ¼ μðϵÞ þ iωðϵÞ; λðϵÞ ¼ μðϵÞ À iωðϵÞ with ωðϵÞ > 0; and Following the standard procedure given in Supplementary Note 6, we finally obtain the Poincaré normal form given in Eq. (9) from which the frequency and amplitude are finally quantified.
Modulating both components. To approximate the amplitude of the first component u 1 from the normal form, we carefully determine the coefficient of the eigenvector corresponding to λ (see Supplementary Note 6 for details). Also, if we want to implement a commensurate variation on the second component u 2 simultaneously, the transformation acted on u 2 should also be considered. Otherwise, as stated in "Results", the amplitude of u 2 could be increased or decreased, whereas that of u 1 is invariant, or vice versa. Our computation [Supplementary Eq. (S17)] shows that the approximation for the amplitude of u 2 includes a coefficient ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Àða 21 þ f 21 Þ=ða 12 þ f 12 Þ p . We finally deduce Eq. (10) [Supplementary Eq. (S18)] to eliminate this influence.
The F-N system and its coordination. The R-D version of the F-N model is given by where all system parameters are positive and ϵ is the bifurcation parameter. With appropriate parameter values, it possesses a positive stationary solution denoted by (V, W) = (V 0 , W 0 ). Before investigation, we first convert the system into the form of Eq. (11). The formulation of the coordinator together with the detailed computations can be found in Supplementary Note 7.
The index for assessing energy consumption. For different biological models, distinct definitions may be introduced to assess the energy consumption. For the diffusive neuronal model considered here, we adopt the classical definition used in the cable model 71 . Accordingly, the energy at time t* induced by the stimulus current I in a diffusive neuronal model is defined as Since we consider the periodic oscillation in this work, in order to make a fair comparison, we take the average energy consumption into account. It is defined as where T is the period of the considered oscillation. To implement our designed coordinator, the linear terms f 11 (V − V 0 ) + f 12 (W − W 0 ) (Supplementary Note 7) can be integrated as an additional stimulus current I c . Our aim is to assess the energy consumption of the coordinator alone rather than the entire system. We therefore refer to an analogous power consumption index used in controlling the activity of the F-N neuron 72 . Consequently, we finally define the index for our purpose as follows where I o is the original constant stimulus current and V o and V c are the membrane potential in the original and controlled system, respectively. This index evaluates the average energy consumption of the designed coordinator. As for the other biological models such as the "cancer network" considered in this work, the energy can be evaluated by the average L 2 -norm of the designed coordinator over one period, which yields the following index where T is the period of oscillation and ∥ ⋅ ∥ is the L 2 -norm induced by the inner product [Supplementary Eq. (S4)]. We introduce L 2 -norm here because it is a generic definition for the assessment of a control policy 73 including the one used in a biological network with protein-protein interactions 63 . Moreover, the biorhythm can be regarded as a cyclic signal for information processing whose energy is also conventionally defined as the integral of L 2 -norm 74 . Therefore, it is appropriate for the "cancer network". Of course, there could be other definition for a different model. It is worth pointing out that the two definitions introduced here can be computed in practice once the signal of a biorhythm (e.g., action-potential or expression level) is measured.
Numerical simulations. All the numerical simulations are performed using PDE solver pdepe in MATLAB R2019a (The Mathworks), and both absolute and relative tolerance are set to 1e−8.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.