 Research
 Open access
 Published:
An oscillating reaction network with an exact closed form solution in the time domain
BMC Bioinformatics volume 24, Article number: 466 (2023)
Abstract
Background
Oscillatory behavior is critical to many life sustaining processes such as cell cycles, circadian rhythms, and notch signaling. Important biological functions depend on the characteristics of these oscillations (hereafter, oscillation characteristics or OCs): frequency (e.g., event timings), amplitude (e.g., signal strength), and phase (e.g., event sequencing). Numerous oscillating reaction networks have been documented or proposed. Some investigators claim that oscillations in reaction networks require nonlinear dynamics in that at least one rate law is a nonlinear function of species concentrations. No one has shown that oscillations can be produced for a reaction network with linear dynamics. Further, no one has obtained closed form solutions for the frequency, amplitude and phase of any oscillating reaction network. Finally, no one has published an algorithm for constructing oscillating reaction networks with desired OCs.
Results
This is a theoretical study that analyzes reaction networks in terms of their representation as systems of ordinary differential equations. Our contributions are: (a) construction of an oscillating, two species reaction network [two species harmonic oscillator (2SHO)] that has no nonlinearity; (b) obtaining closed form formulas that calculate frequency, amplitude, and phase in terms of the parameters of the 2SHO reaction network, something that has not been done for any published oscillating reaction network; and (c) development of an algorithm that parameterizes the 2SHO to achieve desired oscillation, a capability that has not been produced for any published oscillating reaction network.
Conclusions
Our 2SHO demonstrates the feasibility of creating an oscillating reaction network whose dynamics are described by a system of linear differential equations. Because it is a linear system, we can derive closed form expressions for the frequency, amplitude, and phase of oscillations, something that has not been done for other published reaction networks. With these formulas, we can design 2SHO reaction networks to have desired oscillation characteristics. Finally, our sensitivity analysis suggests an approach to constructing a 2SHO for a biochemical system.
Background
Oscillatory behavior is critical to many life sustaining processes. Examples include: cell cycles [1], circadian rhythms [2], notch signaling in the development of the nervous system [3], tissue development [4], gene transcription [5], and efficient signaling [6]. Biological oscillators are also important elements in building applications in synthetic biology [7,8,9].
The characteristics of biological oscillations often have critical biological functions. Frequency is used to control the times at which events are initiated, such as circadian cycles and chromatin modifications [10]. Amplitude controls the strength of signaling [11]. Phase plays a role in the sequencing of processes within the cell cycle [12]. Since biological oscillators typically cause changes in the concentration of chemical species, these oscillators must have a DC offset so that values are nonnegative. Collectively, we refer to frequency, amplitude, phase, and DC offset as oscillation characteristics (OCs).
One way to understand the relationship between an oscillating reaction network and its OCs is to construct a closedform, timedomain solution of the network’s behavior in terms of parameters such as kinetic constants and initial concentrations of chemical species. From these mathematical expressions, we obtain insights such as: (a) if one or more reactions are unnecessary to achieve oscillations; (b) relationships between the kinetic constants of reactions; and (c) how to assign values to network parameters so as to achieve desired OCs.
Many researchers have investigated structural aspects of oscillating reaction networks. These structures include: positive feedback, negative feedback, balancing reaction rates, and ultrasensitivity [13,14,15,16,17]. Others have built biological oscillators [7, 8, 12, 18,19,20,21]. But neither kind of investigation addresses our interest in a closedform, timedomain solution that relates parameters of a reaction network to its OCs.
More relevant to our work are quantitative models of biological oscillators. For the most part, existing models are systems of nonlinear ordinary differential equations [22,23,24,25,26,27]. Typically, these nonlinearities arise because one or more rate law is a nonlinear function of species concentrations (e.g., Michaelis–Menten kinetics). The complexity of these models prohibits the construction of a closedform, timedomain solution.
We are aware of two approaches that circumvent the limitations of nonlinear ODEs. The first uses an empirical approach, system identification, to construct a linear model that approximates the nonlinear system (e.g., [28]). Models constructed in this way provide accurate predictions near the operating point at which system identification is done. The second approach constructs a linear approximation to a nonlinear ODE (e.g., [29]). Typical approximations make assumptions about relative reaction rates and/or magnitudes of species concentrations. In both cases, the construction of a linear models greatly reduces the complexity of the mathematical expressions, and this in turn makes it possible to obtain a closedform, timedomain solution. However, the approximations limit the extent to which the resulting mathematical expressions provide useful interpretations of how the parameters of the reaction network affect OCs.
The present work is a theoretical study in that we do not build a biological oscillator. Instead, we consider reaction networks as modelled by systems of ordinary differential equations. Gonze and Ruoff [30] claim that oscillations in reaction networks require that “the kinetic rate laws of the reaction mechanism must be sufficiently ‘nonlinear’ to destabilize the steady state.” This claim is echoed by [31] as well. We note that the context for these claims may be for biological implementations of oscillators. Our work is a theoretical exploration of the construction of an oscillator for a reaction network whose dynamics can be described by a system of linear ordinary differential equations.
It turns out that such a construction is possible. Indeed, we are able to construct a harmonic oscillator for a two species reaction network. We refer to this network as a two species harmonic oscillator (2SHO). Because the 2SHO has linear system dynamics, its behavior depends on initial conditions (e.g., [32]).
At first glance, this may seem to be a modest undertaking since it is well known that harmonic oscillators can be constructed from a system of linear differential equations. Specifically, harmonic oscillations result when the Jacobian matrix of the linear system has pure imaginary eigenvalues (e.g, [33]). So, why not just populate the Jacobian such that the reaction network has imaginary eigenvalues?
Such an approach is difficult to apply to reaction networks because of constraints on the system of linear differential equations. Foremost, reaction networks must produce concentrations of chemical species that are nonnegative (since negative concentrations are impossible). This turns out to be a nontrivial constraint to satisfy.
A second challenge is that we cannot instantiate a reaction network from an arbitrary Jacobian matrix since the nature of reaction networks imposes constraints on the Jacobian matrix. In particular, a single reaction may affect multiple entries in the matrix. For example, consider a two species reaction network that has a \(2 \times 2\) Jacobian, \({\textbf {A}} = \{ a_{ij} \}\). Suppose there is a reaction \(S_1 \rightarrow S_2\) with the mass action rate law \(kS_1\). This one reaction affects both \(a_{11}\) and \(a_{12}\); it subtracts \(k S_1\) from \(a_{11}\), and it adds the same amount to \(a_{21}\).
A final challenge is that the choice of rate laws must be biologically plausible. That is, we cannot simply invent a rate law that results in a Jacobian with pure imaginary eigenvalues. In our work, a rate law is biologically plausible if it is used in a published model such as BioModels [34].
The contributions of this paper are:

1.
construction of a harmonic oscillator from a two species reaction network (two species harmonic oscillator (2SHO)), a result that demonstrates that an oscillator can be implemented in a reaction network with no nonlinearity;

2.
obtaining closed form formulas that calculate the frequency, amplitude, and phase of the 2SHO in terms of its parameters, something that has not been done for any published oscillating reaction network; and

3.
development of the parameterizeOscillator algorithm that parameterizes the 2SHO to achieve desired oscillation, a capability that has not been produced for any oscillating reaction network.
We use these results to analyze the sensitivity of 2SHO to variations in values of the kinetic constants; typical varaiabilities are between 1 and 20% (e.g., [35, 36]). Our analysis suggests that with 4 to 6 attempts, there is a 90% probability of constructing a 2SHO for a suitable biochemical system.
Methods
Our method is to propose a two species reaction network whose kinetics can be described by a system of linear differential equations. We solve this system (an initial value problem) to obtain a harmonic oscillator with closed form formulas that relate OCs to parameters of the reaction network. Finally, we develop an algorithm for parameterizing the two species harmonic oscillator to achieve desired OCs.
Two species reaction network with linear dynamics
This section develops a biologically feasible reaction network whose kinetics can be described by a system of linear ODEs. Our approach constructs a harmonic oscillator, a model that has been used elsewhere to explain binding and activation between insulin and insulin growth factor receptors [37]. Indeed, harmonic oscillator behaviors are apparent in parameterizations of the LotkaVolterra and Susceptible Infectious Recovered (SIR) models [38], both of which are widely used in biology.
Figure 1 displays our two species harmonic oscillator (2SHO). The construction of 2SHO is driven by the fundamentals of harmonic oscillators in linear systems, as we now describe. The matrix representation of our linear system is:
where \({\textbf {x}} (t) = \{x_n (t)\}\) is an N dimension vector of time varying of species concentrations; \(\dot{{\textbf {x}}} (t)\) is the time derivative of \({\textbf {x}}(t)\); \({\textbf {A}} = \{a_{ij} \}\) is an \(N \times N\) Jacobian matrix of constants; and \({\textbf {u}}\) is an N dimensional vector of constants that are forced inputs. We want to construct a reaction network that has a sustained oscillation. Since this is a linear system, the oscillations will be sinusoids (e.g., [32]).
From the foregoing, we have the following constraints:

C1: Rate laws in the reaction network are a linear function of the concentrations of \(x_n\)(t).

C2: \(x_n(t) \ge 0\) so that the reaction network is biologically feasible.
We simplify the problem by having \(N=2\) since this is sufficient to obtain oscillations. This means that the eigenvalues of \({\textbf {A}}\) must be pure imaginary numbers (e.g., [39]). Let \(\tau = a_{11} + a_{22}\) be the trace of \({\textbf {A}}\), and \(\Delta = a_{11} a_{22}  a_{12} a_{21}\) be the determinant of \({\textbf {A}}\). The eigenvalues are complex conjugates \(\lambda _1, \lambda _2\) such that
Clearly, we obtain pure imaginary eigenvalues only if \(\tau = 0\) and \(\Delta > 0\). Thus, we add the constraints

C3: \(\tau = 0\), where \(\tau\) is the trace of \({\textbf {A}}\).

C4: \(\Delta > 0\), where \(\Delta\) is the determinant of \({\textbf {A}}\).
With these constraints, \(\lambda _n = \pm \theta i,\) where \(\theta = \sqrt{\Delta }.\)
We use \(k_i \ge 0\) to specify kinetics constants, where i indexes the reaction in the network. Let \(S_n\) be a chemical species whose time varying concentration is the state variable \(x_n(t)\) for \(n=\{1,2\}.\) We start by having twoway interactions between the species. Rate laws are specified above the reaction arrow.

\(R_1: S_1 \xrightarrow {k_1 S_1} S_2\)

\(R_2: S_2 \xrightarrow {k_2 S_2} S_1\)
The foregoing reactions have mass action kinetics, which is widely used in models of chemical systems. With just these reactions, \({\textbf {A}}\) is
Clearly, C3 does not hold unless we eliminate \(R_1, R_2\) by having \(k_1 = 0 = k_2\) so that \(\tau = 0.\) Instead, we make \(a_{11}\) positive by adding an autocatalysis reaction, a kind of reaction that arises in many biological oscillators [40].

\(R_3: S_1 \xrightarrow {k_3 S_1} 2 S_1\)
and so \({\textbf {A}}\) becomes
which is positive with appropriate choices of the \(k_i\).
Since \(R_3\) synthesizes \(S_1\), we need a reaction that degrades \(S_1\) in order for the system to be stable. If we use mass action kinetics for this reaction, \(kS_1\), it changes \(a_{11}\) to \(k_1k_2k\), which makes it more difficult to satisfy C3. An alternative is a fixed degradation rate of \(k_4 > 0\), where \(k_4 = u_1\) is the first element of the vector \({\textbf {u}}\) in Eq. (1).

\(R_4\): \(S_1\xrightarrow {k_4}\emptyset\)
Examples of fixed rate degradation reactions in BioModels are: reaction reaction_0 in BIOMD0000000112, reaction ATP_Jerp in BIOMD0000000059, and reaction inhibition_parameter2 in BIOMD0000000224.
We must still address constraint C4, that the determinant is positive. This requires that \(a_{12} a_{21} < 0,\) which mean that one of \(a_{12}, a_{21}\) is positive and the other is negative. We make \(a_{21} < 0\) by degrading \(S_2\) at a rate controlled by \(S_1\). That is,

\(R_5\): \(S_2 \xrightarrow {k_5 S_1}\emptyset\)
Having the degradation of one chemical species controlled by another chemical species is common in BioModels. Some examples are: the fast reaction in model BIOMD0000000108, reaction r10 in BIOMD0000000145, and reaction RuBisCO_5_EOP in BIOMD0000000392. With the addition of reaction \(R_5\), \({\textbf {A}}\) becomes
The final reaction in our network compensates for degrading \(S_2\) by synthesizing this species at a fixed rate. That is,

\(R_6\): \(\emptyset \xrightarrow {k_6} S_2\)
(We do not explicitly cite examples of similar synthesis reactions since they are widely used in BioModels.) From this, we observe that
C3 and C4 constrain the values of the kinetic constants. From C3, we know that
From C4, we know that \((k_3  k_1)(k_2)  k_2(k_1  k_5) > 0,\) or \(k_3 < k_5.\) We define \(k_d >0\) such that
And so, \(k_5 = k_1 + k_2 + k_d\). This gives us
And from this we calculate the determinant of A:
And hence, the frequency \(\theta\) is
Further, \(k_2, k_d > 0\) implies that
Since \(k_3, k_5\) are calculated from other parameters, we refer to them as dependent parameters. \(k_1, k_2, k_4, k_6, k_d, x_1 (0), x_2 (0)\) are the independent parameters.
A brief technical note on Eq. (7). There are actually two solutions for \(\theta\), \(\pm \sqrt{k_2 k_d}\). These solutions result in oscillations at the same frequency but with different phases. Our approach is to treat phase as a separate oscillation characteristic, and so we just use \(\theta\) for the frequency.
We can now satisfy 3 of our 4 constraints. C1 is satisfied since all rate laws are linear in \(x_n (t)\), the time varying concentrations of \(S_n\). C3 is satisfied since \(\tau = k_2 + ( k_2) = 0.\) And, C4 is satisfied by Eq. (8). To address C2, that \(x_n(t) \ge 0\), we must find the time domain solution of the reaction network.
Time domain solution for two species harmonic oscillator
Solving Eq. (1) is an initial value problem, where the initial values of \(S_1\), \(S_2\) are \(x_1(0), x_2(0)\). We proceed as follows: (a) solve the homogeneous equation \({\dot{{\textbf {x}}}}^H = {\textbf {A}} {\textbf {x}}^H(t)\); (b) find a particular solution such that \(\dot{{\textbf {x}}}^P (t) = {\textbf {A}} {\textbf {x}}^P(t) + {\textbf {u}}\); and (c) properly structure the complete solution \({\textbf {x}}(t) = \begin{pmatrix} x_1 (t) \\ x_2(t) \end{pmatrix} = {\textbf {x}}^H (t) + {\textbf {x}}^P (t)\) so that we isolate terms for amplitude, frequency, phase, and DC offset. The derivation is a bit long, and so full details are reported in the Additional file 1.
The solutions for \(x_n(t)\) have the form
where \(\alpha _n\) is the amplitude of oscillation for \(x_n(t)\), \(\theta\) is the frequency in radians, \(\phi _n\) is the phase in radians, and \(\omega _n\) is the DC offset. \(\alpha _n, \theta , \phi _n, \omega _n\) are functions of the \(k_i\) and \(x_n (0)\). Table 1 displays the formulas for the OCs. Because of technical details related to the inverse tangent function, \(\phi _n\) depends on the term \(\pi _n\). (See the Additional file 1 for more on these technical details.) These terms are:
Similarly,
With the symbolic solution in hand, we return to constraint \(C_2\), \(x_n (t) \ge 0.\) This is equivalent to \(\omega _n \ge \alpha _n.\) Unfortunately, the complexity of the OC formulas makes it difficult to solve these formulas to find constraints on parameters that ensure \(C_2\). Indeed, it is frustratingly difficult to find special cases in which C2 is satisfied. For example, consider the situation in which \(k_2\) is large. This results in the following:
where \(T = \sqrt{k_{4}^{2}  2 k_{4} k_{6} + k_{6}^{2} + \theta ^{2} (x_1 (0))^{2} + 2 \theta ^{2} x_1 (0) x_2 (0) + \theta ^{2} (x_2 (0))^{2}}.\) We can satisfy Eq. (10) (i.e., \(x_1 (t) \ge 0\)) by making \(k_6\) large and \(k_4\) small. However, to satisfy Eq. (11) (i.e., \(x_2(t) \ge 0\)), we must do the opposite: make \(k_6\) small and make \(k_4\) large.
Parameterizing the two species harmonic oscillator
Since it is difficult to construct a symbolic expression for the constraint \(x_n (t) \ge 0\), we proceed numerically by developing an algorithm that finds values of the independent parameters that achieve desired OCs for the 2SHO. Desired OCs are indicated by the superscript \(\star\): \(\theta ^{\star }, \alpha ^{\star }, \phi ^{\star }, \omega ^{\star }\).
We proceed by formulating an optimization problem. The objective is to find values of independent parameters that result in the best match with \(x^{\star }(t) = \alpha ^{\star } sin(\theta ^{\star } t + \phi ^{\star }) + \omega ^{\star }\). This is done by searching the space of possible values of independent parameters. We denote elements of the search space by \({\mathcal {P}}\), where \({\mathcal {P}}\) is a vector that specifies the value of each independent parameter. We want to find \(\hat{{\mathcal {P}}}\) in this space such that \(x_n(t; \hat{{\mathcal {P}}})\) is close to \(x^{\star }(t)\) for one of \(n \in \{1, 2\}.\)
Given candidate parameters \({\mathcal {P}}\), we can apply the formulas in Table 1 to calculate \(x_n(t; {\mathcal {P}}).\) The quality of this candidate is assessed using the loss function
This gives us the solutions
To simplify our notation in the sequel, \({\hat{x}}_n (t) = x_n(t; \hat{{\mathcal {P}}_n}).\)
Our algorithm constructs a single sinusoid. Does it matter whether we use \(S_1\) or \(S_2\)? We refer to the chosen species as the chosen oscillating species. It turns out that depending on the desired OCs, sometimes it is better to choose \(S_1\); other times \(S_2\) works better. The algorithm selects the chosen oscillating species based on which loss function is smaller, \(L_1 (\hat{{\mathcal {P}}}_1)\) or \(L_2 (\hat{{\mathcal {P}}}_2)\). We denote the chosen oscillating species by \({\hat{x}} (t).\)
We observe that the optimization can be simplified in a couple of ways. First, the only constraint on \(k_1\) is \(k_1 \ge 0.\) So, we set \(k_1 = 1\). Second, we can calculate \(k_d\) from the desired frequency \(\theta ^{\star }\) using Eq. (7), and so \(k_d = \frac{(\theta ^{\star })^2}{k_2}\). The combination of these simplifications reduces the dimension of the search space from 7 to 5.
Figure 2 displays parameterizeOscillator, our algorithm for designing an oscillator for the reaction network in Fig. 1. The python implementation of the algorithm is in the module designer.py in the github repository for this project. The implementation uses the python package lmfit to find \(\hat{{\mathcal {P}}}_n (\theta ^{\star })\) using gradient descent.
There are some important details in the implementation of parameterizeOscillator. First, the quality of the optimization is greatly improved by doing multiple calls to lmfit using randomly chosen initial values for the independent parameters to more effectively scan the search space. Second, in calculating the loss functions, it is critical to adjust the density of time values (t) to the value of \(\theta ^{\star }\) so that the algorithm sees at least two complete sinusoids at 10 or so different phases. Last, properly handling of the hard constraint \(x_n(t) \ge 0\) is required so that gradient decent works effectively. We use relaxation [41], a technique that treats the hard constraints as soft constraints with large weights. This means that instead of using \(L_n(({\mathcal {P}})\) in Eq. (12) as the loss function, we use \(L^+_n({\mathcal {P}})\):
where \(1_x\) is 1 if \(x < 0\) and is 0 otherwise, and w is a large number. \(L_n^{+} ({\mathcal {P}})\) has large gradients when \(x_1(t) <0\) or \(x_2(t) < 0\) so that gradient descent moves away from regions of the search space in which \(C_2\) is violated.
Observe that parameterizeOscillator could be implemented using simulation instead of the formulas in Table 1. However, doing so would require considerably more computational resources. First, simulating the reaction network is several orders of magnitude slower than evaluating the formulas in Table 1. Second, the simulation approach does not exploit the relationships between parameters that reduce the size of the search space (e.g., calculating \(k_d\) from \(k_2, \theta ^{\star }\)). Thus, a simulationbased algorithm would have a search space consisting of seven parameters in contrast to the 5 parameters used with the formulas based approach. This larger search space requires considerably more computational resources.
Results
This Section analyzes the characteristics of the 2SHO and related analysis developed in the previous section.
Accuracy of predictions of oscillation characteristics
We begin by investigating the accuracy of predictions of 2SHO oscillation characteristics that are obtained from the formulas in Table 1. Figure 3 displays simulations of the 2SHO network detailed in Fig. 1 for four different values of the independent parameters. (Dependent parameters are calculated as described above.) Lines in the plots are simulation results, and the markers are model predictions using the formulas in Table 1. We see that in all cases, model predictions coincide with simulations results. We have done thousands of such simulations. In all cases, model predictions are identical with simulation results (within a tolerance of \(10^{13}\) to account for numerical errors). Although this does not prove the correctness of the model, it is strong confirmation. The ultimate proof is the correctness of our derivations as detailed in the Additional file 1.
Next, we investigate our algorithm for parameterizing the 2SHO as detailed in Fig. 2. Recall that the inputs to the algorithm are desired oscillation characteristics \(\theta ^{\star }, \alpha ^{\star }, \phi ^{\star }, \omega ^{\star }\) that generate the sinusoidal concentrations \(x^{\star }(t) = \alpha ^{\star } sin(\theta ^{\star } t + \phi ^{\star }) + \omega ^{\star }\). The algorithm finds independent parameters \(\hat{{\mathcal {P}}} (\theta ^{\star })\) that generate the concentrations \({\hat{x}} (t)\) for one of \(S_1, S_2\) with the constraint that there are no negative concentrations.
We use the term design error to refer to deviations between \(x^{\star }(t)\) and \({\hat{x}}(t)\). We consider several kinds of design errors. A feasibility design errors occurs if \({\hat{x}}_n (t) < 0\) for some n, t. An amplitude design error is a deviation from \(\alpha ^{\star }\). Let \(\hat{\alpha }\) be the amplitude achieved by \({\hat{x}}(t).\) The amplitude design error is calculated as \(\frac{\hat{\alpha }  \alpha ^{\star }}{\alpha ^{\star }}\). Similarly, phase design error is the deviation from \(\phi ^{\star }\) of the phase (as a fraction of a cycle) achieved by \(\hat{{\mathcal {P}}}(\theta )\). The phase design error is \(\frac{\hat{\phi }  \phi ^{\star }}{2\pi }\). Note that we do not consider errors in the frequency that result from the parameters returned by the design algorithm. This is because the algorithm uses Eqs. (7), (3), and (4) to ensure that there are oscillating concentrations at the desired frequency \(\theta ^{\star }\).
We simplify our studies by setting the desired DC offset \(\omega ^{\star }\) to the desired oscillation amplitude \(\alpha ^{\star }\). Our studies are conducted over 3 decimal orders of magnitude for both frequency (radians/sec) and amplitude: \(\theta ^{\star }, \alpha ^{\star } \in [0.1, 100]\) in 8 increments. We consider four values of phase that are likely the most problematic because of challenges with calculating the inverse of the tangent function: \(\phi ^{\star } \in \{0, \frac{\pi }{2}, \pi , \frac{2 \pi }{3} \}.\) Throughout, we set the maximum value of the kinetic constants (\(k_i\)) to 1000 to make the results more credible, but doing so can increase design errors.
We begin with feasibility design errors. parameterizeOscillator is very robust to feasibility design errors. We have conducted several thousand simulations, and only a couple of them returned values of independent parameters that had negative values for the concentrations of \(S_1, S_2\).
To analyze amplitude design errors, we first consider a variant of the algorithm in Fig. 2 in which the algorithm always returns \(\hat{{\mathcal {P}}}_1 (\theta )\), the parameter estimates obtained if \(S_1\) is the chosen oscillating species. Figure 4 displays the results of these studies. There are four plots, one for each value of phase. Each plot is a heatmap with desired frequency (\(\theta ^{\star }\)) as the horizontal axis and desired amplitude (\(\alpha ^{\star }\)) as the vertical axis. Cells are colored by the magnitude of the design error, and they are annotated with the value of the design error followed by a letter. The letter “a” means that \(S_1\) is the chosen oscillating species; “b” indicates that \(S_2\) is the chosen oscillating species. Only the letter “a” appears in these plots since \(S_1\) is always the chosen oscillating species.
Amplitude design error is in the range \([1, \infty ]\). A 0 means that there is no design error; a 1 means that the chosen oscillating species always has a concentration of zero (and so there is no oscillation). Amplitude design error is mostly 0 in the plots, except when both amplitude and frequency are large. The reason is that in our studies, the maximum value of kinetic constants is 1000. When both \(\theta ^{\star }\) and \(\alpha ^{\star }\) are large, much larger values of \(k_4, k_6\) are required. Figure 5 plots the results of studies in which either \(S_1\) or \(S_2\) may be the chosen oscillating species, as is done in parameterizeOscillator. We see that there is a significant reduction in amplitude design error.
Figure 6 displays the results for phase design error (in units of the fraction of an oscillation cycle). We see that phase design errors are mostly 0, although occasionally there is an error of 0.1 or \(\) 0.1.
Figure 7 studies the distribution of parameter values for the foregoing studies. Note that \(k_1=1\) by design since from Table 1, we know that \(k_1\) does not influence the behavior of the reaction network. Also, by design, the maximum value of a kinetic constant is 1000.
The parameters \(k_2\), \(k_3\), and \(k_5\) are mass action kinetic constants for reactions with a single reactant. We see that they mostly have small values, although there are some instances in which these constants exceed 500. On the other hand, the zeroth order kinetic constants \(k_4\) and \(k_6\) tend to be much larger. Because our studies restrict parameterizeOscillator to choose \(k_4, k_6\) with values less than 1000, there are larger amplitude design errors for studies in which both the frequency and amplitude are large.
Parameter sensitivity and plausibility of biochemical implementation
Here we explore the sensitivity of the oscillation characteristics of 2SHO to variation in parameter values, especially with the thought of a biochemical implementation of the harmonic oscillator.
We start by analyzing the relative sensitivity of frequency, \(\theta\), to variation in values of a parameter k. This relative sensitivity is denote by \(S_k^{\theta }\). As defined in [42]:
Note that this is a unitless quantity. From Eq. (7), we know that \(\theta = \sqrt{ k_2 k_d} = \sqrt{k_2 (k_5  k_3)}\), and so:
\(S^{\theta }_k = 0\) for \(k \notin \{k_2, k_3, k_5 \}\). To put this in perspective, we likely want \(k_5>> k_3\) to ensure that we get oscillations. So, for \(k \in \{k_2, k_3, k_5 \}\), \(S^{\theta }_k \approx \frac{1}{2}\). That is, there is a 0.5% change in \(\theta\) for a 1% change in k.
The remaining OCs are complicated functions of the 2SHO parameters, and so are difficult to analyze analytically. Instead, we conduct numerical studies in which values of the 2SHO parameters are chosen randomly. The primary objective of these studies is to determine the feasibility of constructing a 2SHO for a biochemical system. The details of the biochemistry are beyond the scope of this paper. However, we do address various considerations for selection of the species and reactions since 2SHO imposes some constraints.
The first constraint comes from Eq. (3), that \(k_3 + k_1 + k_2\). This constraint is difficult to achieved by a priori design because of the variability of values of kinetic constants. We propose addressing this constraint by having reaction \(R_3\) be controlled enzymatically and operate in a region of low enzyme concentrations so that mass action kinetics are a good approximation. By titrating the enzyme, it should be possible to be fairly close to the Eq. (3) constraint. We refer to this as the titration approach.
A second constraint imposed by 2SHO is that \(k_3 < k_5\). This constraint is easily realized by our choice of 2SHO parameters. Specifically, we use: \(k_1=1.00,~ k_2=3.91,~ k_3=4.91, k_4=4.77,~ k_5=15.00,~ k_6=92.2,~ x_1(0)=5.00,~ x_2(0)=10.49\). This produces the OCs \(\theta =6.28,~ \alpha _1=3.00,~ \alpha _2=5.67,~ \phi _1=0.00, \phi _2=2.12,~ \omega _1=5.00,~ \omega _2=5.67\). Clearly, \(k_5 \gg k_3\).
The third constraint is that concentrations must be nonnegative; that is, \(\omega _n \ge \alpha _n\). Since this is a nontrivial constraint to assess, we analyze it numerically. We consider the values of the 2SHO parameters to be random variables whose distributions are centered on the parameters in the previous paragraph. Many have explored the experimental error associated with estimating kinetic constants (e.g., [35, 36]). From this literature, it seems that the standard deviations of these errors are in the range of 1–20% of the value of the kinetic constant. Thus, our studies use a normalized standard deviation (also known as the coefficient of variation), which is the measured standard deviation divided by the sample mean.
Figure 8a displays our results for feasibility, those situations when the randomly chosen parameters satisfy the constraints that \(k_5 \ge k_3\) and \(\omega _n > \alpha _n\). (Vertical bars \(are \pm two\) standard deviations.) Feasibility is about 0.57 when normalized standard deviation is 0.01, and it declines to about 0.33 when normalized standard deviation is 0.2. Put differently, the probability p(m, f) of making a 2SHO in m tries with a feasibility of f for each try is \(p(m,f) = 1(1f)^m\). Note that \(p(4; 0.57) \approx 0.9 \approx p(6; 0.37)\). That is, based on the variability of estimates of kinetic constants, there is a 90% chance of making a 2SHO in 4 to 6 attempts.
Figure 8b–h display measures of the deviation of the oscillation characteristics of situations in which a feasible oscillator is obtained. For \(\alpha _n, \omega _n\), we use the absolute value fractional deviation, \({\hat{z}}z/z\), where \({\hat{z}}\) is the randomly chosen value that is centered at z. This doesn’t work for \(\phi\), since \(\phi _1=0\). So, for \(\phi\), we do not divide by the desired value. As with feasibility, we see that the deviation from desired OCs increases with normalized standard deviation. Of the OCs, frequency \(\theta\) seems to be affected the least.
Discussion
This section explores the 2SHO in more detail using the formulas in Table 1.
We begin by examining the parameters of the reaction network. Consider \(k_1\), the kinetic constant for reaction \(R_1\). \(k_1\) does not appear in Table 1. Indeed, the only reference to \(k_1\) is in Eq. (3) to calculate \(k_3\). So, \(k_1\) can be any nonnegative number. If \(k_1= 0\), we have effectively eliminated reaction \(R_1\). We have done many thousands of simulations in which \(k_1 = 0\) with various values of the independent parameters and calculating the dependent parameters as described above. In all cases, we obtain the oscillating networks predicted by Table 1. From this we conclude, that \(R_1\) is not required.
Next, consider \(k_2\), the kinetic constant that controls the rate at which \(S_2\) is converted to \(S_1\). That is, because of reaction \(R_2\), if \(S_2\) is larger at time t, then \(S_1\) will be larger at a future time, say \(t + t_d\) for \(t_d >0\). In essence, \(k_2\) controls the phase shift from \(S_2\) to \(S_1\). We can see this using the results in Table 1. If \(k_2\) is large, then the two species have the same phases: \(lim_{k_2 \rightarrow \infty } \phi _1 = \frac{k_4 k_6}{x_1(0) + x_2(0)} = lim_{k_2 \rightarrow \infty } \phi _2\). This phase shift in turn determines the frequency, since \(\theta = \sqrt{k_2 k_d}.\) So, a large rate at which there is a change in phase in turn results in higher frequency oscillations.
Now consider \(k_3, k_5\), the kinetic constants for reactions \(R_3\) and \(R_5\). \(R_3\) is an autocatalysis reaction in which \(S_1\) produces two copies of itself. This is a kind of positive feedback at the rate \(k_3 S_1\). \(R_5\) is degradation reaction in which \(S_2\) is eliminated at the rate \(k_5 S_1\). This is a kind of negative feedback in that a larger concentration of \(S_1\) reduces the concentration of \(S_2\), which in turn reduces the concentration of \(S_1\) (because of \(R_2\)). From Eq. (4), we know that \(k_5 > k_3\). Indeed, when \(k_5 = k_3\), then the only eigenvalue of the system is 0. So, negative feedback must be larger than positive feedback in our 2SHO.
Reactions \(R_4, R_6\) have zeroth order kinetics. These are essentially external tuning knobs that adjust oscillation characteristics in complex ways. An extreme case is in the DC offsets \(\omega _n\): if \(k_4 = 0 = k_6\), then \(\omega _1 = 0 = \omega _2\). Note that \(k_4, k_6\) appear in every formula for the OCs in Table 1 except frequency \(\theta\).
Finally, we address the initial concentrations of \(S_1, S_2\), which are denoted by \(x_n(0)\). We see that initial concentrations affect amplitude, as is expected for a harmonic oscillator. We see that initial conditions also affect phase. However, initial concentrations do not impact frequency or DC offset.
Next we comment on related work in light of our methods and results. We start with claims related to nonlinearity. Gonze and Ruoff [30] claims that “the kinetic rate laws of the reaction mechanism must be sufficiently ‘nonlinear’ to destabilize the steady state.” This claim is echoed by [31] as well. Our 2SHO is a counter example to these claims in that its kinetics are described by a system of linear differential equations.
We construct 2SHO from two species because this simplifies the analysis for ensuring that the largest eigenvalue is a pure imaginary number. It is an open question if such a construction can be done for three or more species, and if so, can we obtain closed form solutions for the oscillation characteristics.
Another remark worthy of comment is that oscillations can result from a combination of positive and negative feedback [31]. Indeed, Fig. 1 has both positive and negative feedback. Reaction \(R_3\) provides positive feedback at a rate \(k_3 S_1\) through autocatalysis, and reaction \(R_5\) provides negative feedback at the rate \(k_5 S_1\) by degrading \(S_2\). We add to these remarks the observation that the rate of negative feedback must be larger than the rate of positive feedback, at least in our reaction network, since by Eq. (4), \(k_5 > k_3.\)
Our final remark is a bit more speculative. Disciplines such as electrical and mechanical engineering make extensive use of system identification and control theory in their designs and analyses. A key element of these techniques is the use of frequency analysis such as Bode plots. Frequency analysis requires the ability to generate sinusoids with specific oscillation characteristics, especially frequency, amplitude, and phase. While others have demonstrated the construction of biological clocks (e.g., [43, 44]), there is no technique for generating sinusoids with arbitrary oscillation characteristics. Our parameterizeOscillator algorithm provides a way to choose values of the parameters of the reaction network in Fig. 1 to achieve desired The implementation of a biological signal generator is beyond the scope of this paper. However, such an implementation likely requires enzyme engineering techniques such as those discussed in [40].
Conclusions
Oscillatory behavior is critical to many life sustaining processes, such as cell cycles, circadian rhythms, and notch signaling. We use the term oscillation characteristics (OCs) to refer to the frequency, amplitude, phase, and DC offset of oscillations. OCs affect many biological functions. Examples include: the timing of event initiations can be controlled by oscillation frequency; the strength of signaling can be regulated by the amplitude of oscillations; and the sequence of events can be determined by the phase of oscillations.
This paper is a theoretical study that shows that sustained oscillations are possible in a reaction network with linear dynamics. By linear dynamics is meant that the rate laws are linear in species concentrations. It has not been obvious to many that such networks can produce oscillations. Gonze and Ruoff [30] claim that oscillations in reaction networks require that “the kinetic rate laws of the reaction mechanism must be sufficiently ‘nonlinear’ to destabilize the steady state. ” This claim is echoed by [31] as well. We note that the context for these claims may be for biological implementations of oscillators. Even so, the question remains, can oscillations be generated by a reaction network with linear dynamics?
Although it is fairly starightforward to create a system of linear ODEs that oscillates, it is nontrivial to create an oscillating reaction network with linear dynamics. The challenges are: (a) species concentrations must be nonnegative; (b) the semantics of reactions imposes constraints on the Jacobian matrix (e.g., a single reaction can affect multiple entries in the matrix); and (c) reaction rate laws must be biologically plausible.
Our first contribution is the construction of a reaction network that implements a two species harmonic oscillator (2SHO). This is shown in Fig. 1. We are unaware of others who have implemented a harmonic oscillator in a reaction network. The dynamics of 2SHO are described by a system of linear differential equations. We note that [30, 31] claim that oscillating reaction networks require nonlinearities. This may be true for biological implementations, but it certainly is not true in theory for a reaction network.
Our second contribution is obtaining closed form formulas that calculate the oscillation characteristics (OCs) frequency, amplitude, and phase of the 2SHO in terms of its parameters. The formulas are displayed in Table 1. We know of no other oscillating reaction network for which such formulas have been obtained.
Our third contribution is the development of an algorithm that finds values of parameters that result in a reaction network with desired OCs. The algorithm is detailed in Fig. 2. The algorithm makes use of the OC formulas. We show that the algorithm performs well over a wide range of OCs.
We further analyze the sensitivity of 2SHO to variations in values of the kinetic constants; typical varaiabilities are between 1 and 20%. Our analysis suggests that with 4 to 6 attempts, there is a 90% probability of constructing a 2SHO for a suitable biochemical system.
The OC formulas also allow us to comment on and/or extend other studies of oscillatory reaction networks. One theme in the literature is that nonlinear dynamics are required to create oscillations in reaction networks. Our 2SHO is a counter example to these claims in that its kinetics are described by a system of linear differential equations. Another observation in the literature is that an oscillating network can be constructed by combining positive and negative feedback, an approach used in 2SHO as well. Our formulas extend this insight by showing that for 2SHO, the rate of negative feedback must exceed the rate of positive.
We are very interested in exploring how our results for designing an oscillating linear network might be applied to designing an oscillating nonlinear network. One direction here is using 2SHO to construct a linear approximation to the nonlinear network. If we can develop an appropriate mapping between the two networks, we can tune the linear network in a desired way, and then map these adjustments back to the nonlinear network.
Availability of data and materials
Supplementary material is available at https://github.com/ModelEngineering/Oscillators/blob/main/docs/supplemental_material/supplemental_material.pdf, and the github project for this paper is at https://github.com/ModelEngineering/Oscillators.
Abbreviations
 2SHO:

Two species harmonic oscillator
 OC:

Oscillation characteristic
 ODE:

Ordinary differential equation
 \({\textbf {A}}\) :

Jacobian matrix
 \(\alpha _n\) :

Amplitude of oscillation for species n
 \(\alpha _n^{\star }\) :

Desired amplitude
 \(\hat{\alpha }_n\) :

Amplitude of the oscillation produced by parameterizeOscillator
 \({\textbf {c}} = \begin{pmatrix} c_1 \\ c_2 \end{pmatrix}\) :

Constants associated with the homogeneous solution
 \(\Delta\) :

\(det {\textbf {A}})\) (determinant)
 \({\textbf {F}} (t)\) :

Fundamental matrix
 i :

Indexes kinetic constants
 \(k_i\), \(k_d\) :

Nonnegative constants for kinetic constants
 \(\lambda _n\) :

Eigenvalue
 n :

Indexes states (chemical species)
 N :

Number of species
 \(\omega _n\) :

DC offset of species n
 \(\omega _n^{\star }\) :

Desired DC offset
 \(\hat{\omega }_n\) :

DC offset of the oscillation produced by parameterizeOscillator
 \({\mathcal {P}}\) :

The set of parameters and associated values for the reaction network in Fig. 1
 \({\mathcal {P}} (\theta )\) :

The set of parameters with \(k_d\) determined by \(\theta\) and \(k_2\)
 \(\phi _n\) :

Phase in radians for species n
 \(\phi _n^{\star }\) :

Desired phase
 \(\hat{\phi }_n\) :

Phase of the oscillation produced by parameterizeOscillator
 \(S_n\) :

Chemical species in the reaction network
 t :

Time
 \(\tau\) :

\(tr({\textbf {A}})\) (trace)
 \(\theta\) :

Frequency in radians
 \(\theta ^{\star }\) :

Desired frequency
 \(\hat{\theta }\) :

Frequency of the oscillation produced by parameterizeOscillator
 \({\textbf {u}} = \begin{pmatrix} u_1 \\ u_2 \end{pmatrix}\) :

Forced input (kinetic constants for zeroth order rates)
 \(\dot{\textbf {x}} (t)\) :

Derivative w.r.t. time of \({\textbf {x}}\)
 \(x_n(t)\) :

Time varying concentration of species n
 \(x_n^{\star }\) :

Desired signal
 \({\hat{x}}_n\) :

Signal produced the oscillator constructed parameterizeOscillator
References
Murray AW. Cell cycle extracts. Methods Cell Biol. 1991;36:581–605.
Capper S. Overview of circadian rhythms. Alcohol Res Health. 2001;25(2):85.
Wang R, Liu K, Chen L, Aihara K. Neural fate decisions mediated by transactivation and cisinhibition in notch signaling. Bioinformatics. 2011;27(22):3158–65. https://doi.org/10.1093/bioinformatics/btr551. (ISSN 13674811, 13674803).
Goodwin BC, Cohen MH. A phaseshift model for the spatial and temporal organization of developing systems. J Theor Biol. 1969;25(1):49–107. https://doi.org/10.1016/S00225193(69)800172. (ISSN 10958541).
Nelson DE. Oscillations in NF b signaling control the dynamics of gene expression. Science. 2004;306(5696):704–8. https://doi.org/10.1126/science.1099962. (ISBN: 10959203 (Electronic)).
Bae SA, Acevedo A, Androulakis IP. Asymmetry in signal oscillations contributes to efficiency of periodic systems. Crit Rev Biomed Eng. 2016;44(3):193–211. https://doi.org/10.1615/CritRevBiomedEng.2017019658. (ISSN 0278940X).
Elowitz MB, Leibler S. A synthetic oscillatory network of transcriptional regulators. Nature. 2000;403(6767):335–8. https://doi.org/10.1038/35002125. (ISSN 00280836. ISBN: 00280836 (Print) 00280836 (Linking)).
Perry N, Ninfa AJ. Synthetic networks: oscillators and toggle switches for Escherichia coli. Methods Mol Biol. 2012. https://doi.org/10.1007/978161779412417.
Zhang F, Sun Y, Zhang Y, Shen W, Wang S, Ouyang Q, Luo C. Independent control of amplitude and period in a synthetic oscillator circuit with modified repressilator. Commun Biol. 2022;5(1):23. https://doi.org/10.1038/s42003021029871. (ISSN 23993642).
Venkatachalam V, Jambhekar A, Lahav G. Reading oscillatory instructions: how cells achieve timedependent responses to oscillating transcription factors. Curr Opin Cell Biol. 2022;77: 102099. https://doi.org/10.1016/j.ceb.2022.102099. (ISSN 18790410).
Mahrou B, Pirhanov A, Alijanvand MH, Cho YK, Shin YJ. Degradationdriven protein level oscillation in the yeast saccharomyces cerevisiae. Biosystems. 2022;219: 104717. https://doi.org/10.1016/j.biosystems.2022.104717. (ISSN 03032647).
Ball DA, Marchand J, Poulet M, Baumann WT, Chen KC, Tyson JJ, Peccoud J. Oscillatory dynamics of cell cycle proteins in single yeast cells analyzed by imaging cytometry. PLoS ONE. 2011;6:8. https://doi.org/10.1371/journal.pone.0026272. (ISSN 19326203).
Cao Y, Lopatkin A, You L. Elements of biological oscillations in time and space. Nat Struct Mol Biol. 2016;23(12):1030–4. https://doi.org/10.1038/nsmb.3320. (ISSN 15459985).
Hoon HS, Ferrell JE Jr. Ultrasensitivity part III: cascades, bistable switches, and oscillators. Trends Biochem Sci. 2012;100(2):130–4. https://doi.org/10.1016/j.tibs.2014.10.002.Ultrasensitivity. (ISSN 15378276).
Lenz P, SøgaardAndersen L. Temporal and spatial oscillations in bacteria. Nat Rev Microbiol. 2011;9(8):565–77. https://doi.org/10.1038/nrmicro2612. (ISSN 17401526, 17401534).
Li Z, Yang Q. Systems and synthetic biology approaches in understanding biological oscillators. Quant Biol. 2018;6(1):1–14. https://doi.org/10.1007/s4048401701207. (ISSN 20954697).
Tatka LT, Luk W, Elston TC, Hellerstein JL, Sauro HM. Cesium: a public database of evolved oscillatory reaction networks. Biosystems. 2023;224: 104836. https://doi.org/10.1016/j.biosystems.2023.104836. (ISSN 03032647).
Atkinson Mariette R, Savageau Michael A, Myers Jesse T, Ninfa AJ. Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli. Cell. 2003;113(5):597–607. https://doi.org/10.1016/S00928674(03)003465. (ISSN 00928674).
Nakajima M, Imai K, Ito H, Nishiwaki T, Murayama Y, Iwasaki H, Oyama T, Kondo T. Reconstitution of circadian oscillation of cyanobacterial KaiC phosphorylation in vitro. Science. 2005;308(5720):414–5. https://doi.org/10.1126/science.1108451. (ISSN 00368075).
Rosier BJ, De Greef TF. How to make an oscillator. eLife. 2015;4: e12260. https://doi.org/10.7554/eLife.12260. (ISSN 2050084X).
Weitz M, Kim J, Kapsner K, Winfree E, Franco E, Simmel FC. Diversity in the dynamical behaviour of a compartmentalized programmable biochemical oscillator. Nat Chem. 2014;6(4):295–302. https://doi.org/10.1038/nchem.1869. (ISSN 17554349).
Goldbeter A. A minimal cascade model for the mitotic oscillator involving cyclin and cdc2 kinase. Proc Natl Acad Sci. 1991;88(20):9107–11. https://doi.org/10.1073/pnas.88.20.9107. (ISSN 00278424).
Goodwin BC. Oscillatory behavior in enzymatic control processes. Adv Enzyme Regul. 1965. https://doi.org/10.1016/00652571(65)900671. (ISSN 00652571).
Heinrich R, Rapoport SM (1977) METABOLIC REGULATION AND MODELS MATHEMATICAL i . INTRODUCTION since the identification of the fundamental pathways of intermediate metabolism some 30 years ago it has been recognized that the stage has been set for the study of the regulation of metabolis. Prog Biophys Mol Biol 32.
Huxley H. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952. https://doi.org/10.1109/ICCCT2.2017.7972284.
Rumbell T, Kozloski J. Dimensions of control for subthreshold oscillations and spontaneous firing in dopamine neurons. PLoS Comput Biol. 2019;15(9): e1007375. https://doi.org/10.1371/journal.pcbi.1007375. (ISSN 15537358).
Sadeghpour M, VelizCuba A, Orosz G, Josić K, Bennett MR. Bistability and oscillations in corepressive synthetic microbial consortia. Quant Biol. 2017;5(1):55–66. https://doi.org/10.1007/s404840170100y. (ISSN 20954689, 20954697).
Mahrou B, Pirhanov A, Alijanvand MH, Cho YK, Shin YJ. Degradationdriven protein level oscillation in the yeast saccharomyces cerevisiae. BioSystems. 2022;219: 104717 (ISSN 18728324).
Kut C, Golkhou V, Bader JS. Analytical approximations for the amplitude and period of a relaxation oscillator. BMC Syst Biol. 2009;3(1):6 (ISSN 17520509).
Gonze D, Ruoff P. The goodwin oscillator and its legacy. Acta Biotheor. 2021;69(4):857–74. https://doi.org/10.1007/s10441020093798. (ISSN 15728358).
Novák B, Tyson JJ. Design principles of biochemical oscillators. Nat Rev Mol Cell Biol. 2008;9(12):981–91. https://doi.org/10.1038/nrm2530. (ISSN 14710072).
Ogata K. Modern control engineering. 4th ed. Hoboken: Prentice Hall; 2001.
Alon U. An introduction to systems biology: design principles of biological circuits. Boca Raton: Chapman and Hall/CRC; 2006. https://doi.org/10.1201/9781420011432.
MalikSheriff RS, Glont M, Nguyen TVN, Tiwari K, Roberts MG, Xavier A, Vu MT, Men J, Maire M, Kananathan S, Fairbanks EL, Meyer JP, Arankalle C, Varusai TM, Vincent KSL, Li Corina DR, Dass G, Keating SM, Park YM, Buso N, Rodriguez N, Hucka M, Hermjakob H. BioModels15 years of sharing computational models in life science. Nucleic Acids Res. 2020;48:D407–15. https://doi.org/10.1093/nar/gkz1055. (ISSN 13624962).
Cvetanovic RJ, Singleton DL, Paraskevopoulos G. Evaluations of the mean values and standard errors of rate constants and their temperature coefficients. J Phys Chem. 1979;83(1):50–60.
Götz M, Barth A, Bohr SSR, Börner R, Chen J, Cordes T, Erie DA, Gebhardt C, Hadzic MC, A. S, Hamilton GL, Hatzakis NS, Hugel T, Kisley L, Lamb Don C, De Lannoy C, Mahn C, Dunukara D, De Ridder D, Sanabria H, Schimpf J, Seidel CAM, Sigel RKO, Sletfjerding MB, Thomsen J, Vollmar L, Wanninger S, Weninger KR, Pengning X, Schmid S,. A blind benchmark of analysis tools to infer kinetic rate constants from singlemolecule FRET trajectories. Nat Commun. 2022. https://doi.org/10.1038/s41467022330233. (ISSN 20411723).
Kiselyov VV, Versteyhe S, Gauguin L, De Meyts P. Harmonic oscillator model of the insulin and IGF1 receptors’ allosteric binding and activation. Mol Syst Biol. 2009;5(1):243. https://doi.org/10.1038/msb.2008.78. (ISSN 17444292, 17444292).
Koshkin S, Meyers I. Harmonic oscillators of mathematical biology: many faces of a predator–prey model. Math Mag. 2022;95(3):172–87. https://doi.org/10.1080/0025570X.2022.2055424. (ISSN 0025570X, 19300980).
Vladimir VI. Ordinary differential equation. Berlin: Springer; 2006.
Novichkov AI, Hanopolskyi AI, Miao X, Shimon LJW, DiskinPosner Y, Semenov SN. Autocatalytic and oscillatory reaction networks that form guanidines and products of their cyclization. Nat Commun. 2021;12(1):2994. https://doi.org/10.1038/s41467021232069. (ISSN 20411723).
Boyd S, Vandenberghe L. Convex optimization. New York: Cambridge University Press; 2004. https://doi.org/10.1017/CBO9780511804441.
Sauro HM. Systems biology: introduction to metabolic control analysis. Seattle: Ambrosius Publishing; 2019.
Chavan AG, Swan JA, Heisler J, Sancar C, Ernst DC, Fang M, Palacios JG, Spangler RK, Bagshaw CR, Tripathi S, Crosby P, Golden SS, Partch CL, LiWang A. Reconstitution of an intact clock reveals mechanisms of circadian timekeeping. Science. 2021;374(6564):eabd4453. https://doi.org/10.1126/science.abd4453. (ISSN 00368075, 10959203).
Chuang CH, Lin CL. Synthesizing genetic sequential logic circuit with clock pulse generator. BMC Syst Biol. 2014;8(1):63. https://doi.org/10.1186/17520509863. (ISSN 17520509).
Acknowledgements
JLH thanks Herbert S. Sauro for his thoughtful comments, Veronica Porubsky for comments on this draft, and the reviewers for their helpful suggestions (especially the inclusion of the sensitivity analysis)
Funding
This work was supported by the Washington Research Foundation and by a Data Science Environments project award from the Gordon and Betty Moore Foundation (Award #20131029) and the Alfred P. Sloan Foundation (Award #3835) to the University of Washington eScience Institute.
Author information
Authors and Affiliations
Contributions
JLH conceived of the ideas in this paper, did all derivations and analyses, and wrote the paper.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The author declares that he has no competing interests.
Additional information
Publisher's note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1
. Supplementary material: Detailed derivations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Hellerstein, J. An oscillating reaction network with an exact closed form solution in the time domain. BMC Bioinformatics 24, 466 (2023). https://doi.org/10.1186/s1285902305600w
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1285902305600w