Describing function-based approximations of biomolecular systems.

Mathematical methods provide useful framework for the analysis and design of complex systems. In newer contexts such as biology, however, there is a need to both adapt existing methods as well as to develop new ones. Using a combination of analytical and computational approaches, the authors adapt and develop the method of describing functions to represent the input-output responses of biomolecular signalling systems. They approximate representative systems exhibiting various saturating and hysteretic dynamics in a way that is better than the standard linearisation. Furthermore, they develop analytical upper bounds for the computational error estimates. Finally, they use these error estimates to augment the limit cycle analysis with a simple and quick way to bound the predicted oscillation amplitude. These results provide system approximations that can add more insight into the local behaviour of these systems than standard linearisation, compute responses to other periodic inputs and to analyse limit cycles.

Introduction

Design of biomolecular systems can enable applications in agriculture, medicine and manufacturing [1]. Complementarily, analysing how naturally occurring biomolecular interactions determine cellular behaviour is a fundamental problem in biology [2]. Mathematical frameworks are useful for both these objectives. These provide system representations to test and compare different design choices as a guide to the actual implementation. These also help to develop useful insight into how system interactions can combine to generate the overall behaviour. Mathematical models used in these cases are typically complex, both due to their large scale as well as the inherent nonlinearity, making their analysis challenging. Therefore, there is a need to adapt existing mathematical methods as well as develop these and new ones for the study of such problems.

Biomolecular systems are frequently represented as ordinary differential equations, formulated based on the principles of mass action. The variables in these equations are the concentrations of various biomolecular species that evolve in time depending on their interactions with each other. One approach to study these mathematical models is exhaustive numerical computations, which can catalogue all possible system behaviours, complemented with simpler calculations to understand the key underlying principles [3, 4]. Another approach is theoretical, exploiting the inherent structure to infer system behaviour. An example is the theory of monotone systems [5]. Intermediary approaches may also exist using various approximations to understand system behaviour. Indeed, methods of linear systems theory have often been used to characterise biomolecular system behaviour such as impulse responses in bacterial chemotaxis [6] as well as frequency responses in the osmolarity pathway of the yeast Saccharomyces cerevisiae mediated by the mitogen‐activated protein (MAP) kinase cascade [7, 8, 9] and in the galactose metabolic pathway, also in the yeast S. cerevisiae [10]. Another example of such an intermediary approach is the describing function technique [11, 12], where the frequency response of a nonlinear system to a sinusoidal input of a particular amplitude is approximated by the first harmonic of the resulting output response. This is widely used in classical control engineering to estimate limit cycle behaviour as well as to replace nonlinear input–output responses with corresponding linear approximation. These linear approximations can then be analysed using the well‐developed tools of linear systems theory. In fact, this technique has been applied to analyse biomolecular oscillations [13] and to approximate input–output maps in biomedical contexts [14]. These results present important early work in using this technique for biological systems.

There are at least three striking aspects related to using a describing function‐based linearisation to approximate the input–output response of a biomolecular signalling system. One, describing functions naturally allow the analysis of finite amplitude inputs, as opposed to infinitesimal amplitude inputs in the standard linearisation. These may be more relevant to actual biomolecular contexts such as in the experimental studies of frequency response mentioned above [7, 8, 9, 10] and provide additional insight into the system behaviour. Two, nonlinearities typically analysed using describing functions in the classical contexts are static nonlinearities such as saturation or hysteresis. Contrastingly, the same nonlinearities in biomolecular contexts can have a dynamic character because of the underlying nonlinear biomolecular interactions are embedded in the overall system dynamics. Three, there may be error involved in the approximation that may depend both on input frequency as well as inherent system parameters, which may be important to quantify. Given these, the describing function‐based approximation of these input–output responses including the dependence on system parameters and the nature of approximation error is generally unclear.

Here, we aim to approximate these systems and estimate the resultant error. For this, we used the technique of describing functions, analytically, where possible, as well as computationally. We computed approximations for representative systems with input–output responses exhibiting saturation dynamics with different slopes as well as hysteretic dynamics (formed the basis of preliminary investigation [15]). Next, we computed the approximation error and developed a theoretical error bound for these kinds of systems. Finally, we used these error estimates to augment the classical describing function‐based limit cycle analysis by providing a simple way to estimate the range of oscillation amplitudes. These results adapt the existing method of describing functions for the study of biomolecular systems and should be useful both in analysis and design.

Computation of describing function approximation

Describing function method

Consider a nonlinear system governed by equations where is a vector of species concentrations, u is input which can be a system parameter, is a row vector and y is the output and typically is a nonlinear map. We aim to approximate the input–output system from u to y using describing function [11]. To compute the approximation, we set the input to be , where b is the forcing amplitude, is the forcing frequency and is the input bias. The describing function approximation can be obtained by taking the first harmonic of the output in response to sinusoidal input. This is defined as The magnitude and phase of the approximation can be calculated as respectively.

Next, we use this technique to approximate canonical biomolecular signalling systems. Here, (1) is an ordinary differential equation model of the system under consideration obtained using the law of mass action. As, biomolecular signals cannot be negative, we use the constraint , so that the input is always positive.

Example 1: biomolecular system with saturating input–output map

One of the simplest signal transduction mechanisms that is representative of diverse signalling contexts is of a biomolecular species that can exist in two states [3, 16, 17]. These can interconvert among each other (Fig. 1 inset) depending on the input level. Typically, one of these states has biological activity and can serve as the output. Consider a simple model of this with a biomolecular species (A) that can interconvert between two forms ( and ) at certain forward and reverse rates ( and ). In this two‐state model, the input such as temperature or pheromone levels can be modelled as modulating the rate and the output as the concentration of , the active form of the protein A. A mathematical model in the context of biomolecular signal transduction can be obtained using mass action kinetics

Fig. 1

Approximation of the simple biomolecular system and error estimate

Black line is the steady‐state input–output response for parameters , , and . Blue dot is the operating point used for computation of linearisation. Schematic representation of biomolecular system is shown in the inset, Red triangle shaped, blue pentagon shaped and green circle shaped markers are the magnitude plots of analytical, computational describing function‐based approximations and direct linearisation, respectively, for same parameter set, frequency is varied logarithmically in the range 0.1–, Phase plots corresponding to b, Corresponding error plots

Here, the total concentration of the protein is , which is constant in the simplest scenario considered here. Owing to the presence of the term , where the input term and the output multiply each other, this is a nonlinear equation. The input–output response at steady state can be obtained by setting (Fig. 1).

To compute the approximation of the overall response, we set the input to be and describing function approximation is computed using the output obtained numerically from MATLAB ode23s solver in response to the above input from (2). The results of this computation are shown in Figs. 1.

To obtain an analytical approximation, we set [11] in (4) and collect the like terms to obtain expressions for , and (the detailed steps are in supplementary material S1) The describing function approximation is (Figs. 1) The analytical and computational results match well converging to the linearised frequency response (Figs. 1) as with the operating point .

Parametric dependence of approximation: We note that the phase of the approximation is independent of the forcing amplitude b and coincides with that of linearisation. At frequencies lower than crossover frequency , the output is in phase, whereas at higher frequencies the output lags the input by . The magnitude, on the other hand, depends on forcing amplitude b, in an increasing manner (5). This is a counterintuitive result as steady‐state solution of (4) is increasing and then saturating function of . However, itself decreases as b is increased. This is because decreases as the forcing amplitude b increases (5). Increasing either or increases the crossover frequency , making the output phase similar to input phase. In the following two limits: magnitude decreases with . However, in the high‐frequency limit, magnitude increases with , whereas in the low‐frequency limit, the magnitude first increases and then decreases as is increased.

Finally, , the total protein concentration scales the magnitude but does not change the phase.

This presents an overview of how the describing function approximation depends on system parameters. We note that these relations are more accurate than when obtained by standard linearisation and cannot be achieved through numerical simulation.

Approximation error: It is important to quantify the error in any approximation as a measure of its accuracy and to compare different approximation methods. We computed the relative mean‐square error where y(t) is the numerical solution of (4) with sinusoidal forcing and which can be found analytically from (5) or can be computed using (2) through computer simulation. A comparison of relative error for different approximations (Fig. 1) shows that the computational describing function‐based approximation matches reasonably well with the analytical approximation especially at lower frequencies and is better than the approximation obtained by the standard linearisation.

These results provide a describing function‐based approximation of a simple system as well as a quantification of the error involved.

Example 2: biomolecular covalent modification system with switch‐like input–output map

As a second example of a biomolecular signalling system, we investigated a covalent modification scheme that is present in multiple cellular pathways [18]. This is similar to the example considered above, with enzyme‐catalysed interconversion reactions. Depending on parameter regime of operation, the steady‐state input–output response may have different sensitivities. The reaction scheme for this system is The forward and backward conversions between the two forms A and are catalysed by enzymes and , respectively. As such, these reactions have intermediates and , giving overall conservation relations as where is the total substrate concentration and and are the total enzyme concentrations. Using these the mathematical model can be obtained as As before, we model as the input and the concentration of as the output. The input–output response at steady state can have different sensitivities depending on the parameter regime of operation. Two such choices with a low sensitivity and a high sensitivity are shown in Figs. 2. The nonlinearities are distributed through the system and arise due to the principles of mass action. We aim to approximate the overall input–output response with a describing function‐based linearisation. Therefore, while this is similar to a saturation nonlinearity [12], there are inherent dynamics.

Fig. 2

Approximation of the biomolecular covalent modification system and error estimates

Solid‐black line represents the steady‐state normalised input–output map for parameters , , , , , and , corresponding to a less steep transition. Dashed‐black line represents the steady‐state normalised input–output map which corresponds to a steeper transition, with and other parameters are kept constant. Blue dot represents operating point used for computation of linearisation, Blue‐solid line and red‐solid line represent the magnitude plots of computational describing function‐based approximation and linearisation approximation of less steep transition, respectively, frequency is varied logarithmically in the range 0.1–. Inset shows corresponding phase plots, Blue‐ and red‐dashed lines represent magnitude and phase plots for both approximations in steeper transition case, Corresponding error involved in both approximations for both sensitivity regimes

Calculation of approximation: The describing function approximation is obtained computationally (Fig. 2 – blue‐solid line and Fig. 2 – blue dashed line), as described previously. We also computed, for comparison, the standard linearisation (Fig. 2 – red‐solid line and Fig. 2 – red‐dashed line).

Parameter dependence: The describing function‐based approximations in both regimes have a similar low‐pass nature. The DC gain for the higher sensitivity regime (Fig. 2) is higher than that of the lower sensitivity regime (Fig. 2), consistent with the slope of the operating point. For the low sensitivity regime, the describing function‐based approximation is similar to that of the standard linearisation. Interestingly, however, the DC gain in the higher sensitivity regime is lower than that of the corresponding standard linearisation, while the bandwidth is larger. This shows that, for finite inputs, the gain may not be as high as the standard linearisation predicts. Similarly, the bandwidth may not decrease to a large extent. Therefore, analysis of the describing function‐based approximation provides insight into how the system may behave in case of finite inputs. We emphasise that this realisation is not possible through numerical simulations and the describing function technique provides a better approximation than standard linearisation.

Approximation error: On the basis of the computation (Fig. 2), it can be observed that the error in the describing function‐based approximation is lower than that obtained from the standard linearisation in both parameter regimes considered.

Example 3: biomolecular system with hysteretic input–output map

Another example of a biomolecular input–output response, often encountered in developmental contexts, is hysteresis. One of the simplest ways through which this can be achieved is, through the addition of transcriptional positive feedback in the circuit showed in Example 1. The resulting all‐or‐none behaviour has been experimentally observed during maturation of Xenopus oocytes [19, 20]. Here, we consider a simple system to illustrate the hysteretic input–output map (Fig. 3 inset) where the feedback is and there is the conservation law . The addition of the nonlinear ultrasensitive (n > 1) positive feedback can generate a hysteretic response (Fig. 3). In the hysteretic regime, there are three steady states, two of which are stable and one is unstable. As input parameter is varied, a pair of stable–unstable steady states coalesce causing a monostable high or low steady state. Our aim is to develop an approximation of the entire input–output map from to . While the hysteretic map looks like the classical ones analysed, the difference here is that there are inherent dynamics in the hysteresis nonlinearity.

Fig. 3

Describing function approximation and error estimate for hysteretic system

Solid‐red line represents the steady‐state normalised input–output map for parameters , , n = 2, , K = 1 nM, , corresponding to a hysteretic response. Solid‐blue line represents the steady‐state normalised input–output map when with other parameters same as above. Solid‐black vertical line corresponds to the upper limit of the area where describing function approximation is done for higher‐input amplitude . Dashed‐black vertical lines correspond to the upper and lower limits of the area where describing function approximation is done for lower‐input amplitude, , b = 300. Schematic representation of the system is shown in inset, Red‐solid line and dashed line represent the magnitude plot of computational describing function‐based approximation for hysteretic response with higher and lower‐input amplitudes, respectively, and inset shows phase plot, frequency is varied logarithmically in the range –, Blue‐solid and dashed line show the magnitude and phase plots for non‐hysteretic regime with higher and lower‐input amplitudes, respectively, Corresponding error associated with the approximation in both regimes for both higher and lower‐input amplitudes

Calculation of approximation: The describing function approximation is computed by setting the input, , as a biased sinusoid and taking the first harmonic approximation of the output . We analysed parameter regimes where hysteretic response is present as well as the one in which it is not present (Figs. 3). In fact, due to the presence of multiple steady states in the hysteretic regime, it is unclear which point to take for the standard linearisation. For the describing function‐based approximation, however, there is no such ambiguity.

Parameter dependence: The frequency response shows low‐pass filter nature in both regimes. We note that the frequency response in the hysteretic regime can have a low‐frequency part, where the magnitude is constant and the phase is decreasing (Fig. 3 – solid‐red line). This is exactly the characteristic associated with a pure delay (transfer function ). Therefore, this method of approximating a positive feedback loop adds insight to system analysis by providing a direct map to a pure delay as well as an estimate of the quantitative parameters associated with it. This is consistent with other studies relating positive feedback with a delayed action, for example [21]. With lower‐input amplitude, we still observe two humps in the magnitude and phase plot (Fig. 3), possibly owing to the presence of two time scales – that of protein production and covalent modification. An attenuated version of this trend is also visible in the case when hysteretic response is absent (Fig. 3). Although, the higher‐ and lower‐input amplitude responses do not have much difference when there hysteresis is not expected (Fig. 3).

Approximation error: We calculated the error involved (Fig. 3) in approximation for the hysteretic input–output map for different regimes and different input amplitudes chosen in the input–output plot of Fig. 3. The error is higher for higher‐input amplitude in the hysteretic regime owing to the presence of strong nonlinearity.

Analytical error bound

As with any approximation method, it is important to estimate the associated error. As the describing function technique essentially approximates a dynamical system by the first harmonic of its output, when a sinusoidal input is applied, we develop an error bound using methods associated with Fourier series. The Fourier series of a periodic signal is This can also be represented in exponential form where From the convergence of Fourier series, we know that signals that are continuous or have a finite and bounded discontinuity over a period satisfy [22] where is the Fourier series taken up to harmonic. For the describing function case, we need to estimate which is basically a scaled version of the error defined in (6). For Fourier series, an upper bound of this error can be obtained using the principle of bounded variation [23].

[24]: Let be a function and be a partition of [a, b]. We define the total variation of f over the [a, b] as where the supremum is taken over all partitions of f. The function is said to have bounded variation if is finite in [a, b] and we write .

The following properties are useful in calculating the total variation [24]:

If is monotone in [a, b], then and .

Let be of bounded variation in [a, b]. Then for any and .

Let and c is an arbitrary point in [a, b]. Then if and only if and . Furthermore, .

[23]: If f(t) is periodic with frequency and the total variation over one period is bounded by V, then the mean‐square error in approximating up to harmonic is We use this result to develop an error bound for the describing function approximation. We apply this to classically known static nonlinearities and then illustrate how this can be applied to the dynamic nonlinearities such as the biomolecular systems discussed in Section 2.

Example 4: static saturation nonlinearity

Consider a saturation nonlinearity with a and k denoting the range and slope of saturation. To obtain the describing function approximation, we use a sinusoidal input to the nonlinearity and calculate the first harmonic [12] of the output. The output of the nonlinearity is described as Using the properties of bounded variation In fact, for any real periodic signal, the total variation over a period is bounded and is four times of the amplitude of the signal. Therefore for approximation up to first harmonic. We find that the calculated mean‐squared error in the first harmonic approximation is bounded by the theoretical error bound for static saturation nonlinearity as frequency is varied (Fig. 4).

Fig. 4

Error in describing function approximation and corresponding error bounds

. Blue‐solid line and red dotted line represents mean‐squared error in describing function approximation and corresponding error bound, respectively, for static saturation, when is varied. The same is shown in decibels scale in the inset, Similar plots for static hysteresis case, For the input–output map in Example 1, solid‐blue line represents the mean‐square error in the describing function approximation. Red‐dashed line represents the theoretical error bound when total variation is taken as . Dashed‐blue line represents the bound when total variation is calculated computationally

Example 5: static hysteresis nonlinearity

A relay exhibiting hysteresis is another commonly encountered static nonlinearity [12]. Parameters and D denote the hysteretic angle and the range of saturation, respectively. It works like a simple relay with a phase shift of . The error bound in this case is . Fig. 4 shows the mean‐square error between actual output and the first harmonic approximation, which is bounded by theoretical error bound for all . It can be noted that the error bound for both cases are similar if the maximum magnitude of the output is set to be same in both cases, i.e. ka = D. However, it is observed that the maximum mean‐squared error for the relay with hysteresis is more than that of static saturation nonlinearity. The maximum mean‐squared error increases for saturation case if the slope of saturation increases and tends to the error of relay with hysteresis when k approaches .

Example 1: continued

The nonlinearities encountered in biomolecular systems often have a dynamic character, as analysed in the previous section. To find an error bound for these kinds of nonlinearities, we first need to determine the total variation. We illustrate this for the first example, where we bound the maximum variation for the output to be , the total concentration. Furthermore, the steady‐state input–output response for static saturation case does not take negative values. Comparing these, we find that the total variation of periodic output in this case is and the error bound is . Computationally, we find the total variation from the time trajectory of simulated output. The theoretical as well as computational error bounds for describing function approximation are shown in Fig. 4. Both of these bound the mean‐square error of describing function approximation.

These results provide an error bound for describing function approximation using the concept of total variation and a general way, using the maximum allowable concentration of the output species, obtain these error bounds for biomolecular systems.

Error estimates in limit cycle analysis

Finally, we investigate the use of the error analysis developed above in the analysis of limit cycles, a classically important application of the describing function technique. Although this technique is widely used in control engineering for limit cycle prediction, a treatment of errors involved is relatively less common [25]. Here, we aim to incorporate the simple error estimate developed above to augment the standard analysis.

Example 6: Van der Pol oscillator

We begin with the Van der Pol oscillation, a benchmark nonlinear oscillator and commonly used to illustrate the describing function‐based limit cycle analysis [12]. This oscillator was first analysed in the context of vacuum tube triode circuits and termed as ‘relaxation oscillation’ by Van der Pol [26]. This class of oscillators has multiple applications in physical as well as biological sciences such as in the Fitzhugh–Nagumo oscillator [27, 28]. The dynamical equations of Van der Pol oscillator are This can be separated in a linear and a nonlinear part as shown in Fig. 5.

Fig. 5

Limit cycle analysis of Van der Pol oscillator using describing function technique

Block diagram of Van der Pol oscillator. Solid‐blue line represents nominal describing function approximation, whereas upper and lower limits for the approximation is represented by cyan and red lines, respectively, Black line represents the actual solution for Van der Pol oscillator. Blue, red and cyan line correspond to the solution considering nominal, upper and lower limits of describing function approximation for , Same is repeated for

Describing function approximation with error estimate: We computed the describing function of the nonlinear part (Fig. 5). This is denoted as N. We also computed the error (e) in this approximation to get an uncertainty model , where the  +  and – signs give, respectively, the upper and lower limits. Next, we used the condition for sustained oscillation, 1 + NG = 0, where G is the transfer function of the linear part of the loop. We did this for all three cases – nominal (N), upper limit , lower limit – to get the corresponding values of the oscillation amplitude (A). The Nyquist plot corresponding to the nominal describing function approximation touches the point at rad/s for A = 2. This is consistent with the classical describing function analysis that predicts an oscillation frequency of and an oscillation amplitude of 2. Incorporating the error information in the condition of oscillation, we obtain two other values of oscillation amplitude, A = 1.529 and 3.718. These suggest that the actual oscillations might lie in this uncertainty band.

To check this, we computed the actual solution of the Van der Pol oscillator along with the solution obtained from the nominal as well as the error analysis for and (Figs., 5c and d). At lower values of , the actual solution matches the nominal describing function prediction, though the extent of deviation increases as increases. In fact, the upper and lower amplitudes provide bounds for the variation of the actual amplitude over a period for these parameter choices. Therefore, these error estimates can be used in predicting an uncertainty band of oscillation amplitudes, providing a larger range of parameters where the describing function approximation can help in the limit cycle analysis.

We note that [25] has also performed a general error calculation of limit cycle, which in this case yields error bounds on amplitude and frequency for a larger range of parameter values. Their approach is based on analysing approximation error and its propagation in the entire closed loop. Here, we focus on the bounds simply obtained due to the uncertainty representation of the nonlinearity.

Example 7: biomolecular ring oscillator

The repressilator is a biomolecular oscillator consisting of three proteins that repress each other in a ring [29]. It is a benchmark oscillator in the field of synthetic biology – design using biomolecular substrates – for its demonstration of engineering rich dynamic behaviour inside cells. The dynamical equations for the repressilator are The form of these equations and the nominal parameters are taken from [29]. Briefly, and (i = 1, 2, 3) represent the concentrations of mRNA and protein, respectively. For simplicity, the production and degradation parameters are assumed to be same for all three mRNA–protein pairs. and represent the degradation rate constants for protein and mRNA, respectively. is the translation rate constant. The repression function for the mRNA (), models the repression of as increases (j = 2 for i = 1, j = 3 for i = 2 and j = 1 for i = 3). This is a special case of a general cyclic gene regulatory network which may give rise to oscillations and the existence criteria and general oscillation profiles are given in [30] using harmonic balance analysis. To apply the describing function technique and associated error analysis, we separate the system to linear and nonlinear parts (Fig. 6). There are three transcriptional modules, each of which has one linear and one nonlinear part which is (j = 1, 2, 3). The gain of the describing function approximation for a single nonlinear element is denoted as N. The input to each nonlinear part is concentration of a protein and thus it is set to a biased sinusoidal input, . Unlike the previous case, the input has a bias part. This is because the concentrations in biomolecular systems concentrations cannot be negative .

Fig. 6

Limit cycle analysis of repressilator using describing function technique

Block diagram of repressilator. Solid‐blue line represents nominal describing function approximation, whereas upper and lower limits for the approximation is represented by solid‐red and dashed‐red line, respectively, for parameter values , , , , , K = 40 nM and n = 2, Black, blue and red lines represent the actual solution, nominal describing function solution and describing function solution considering the upper error limit for repressilator when and other parameters are same as b, Same is repeated for

Calculation of frequency and describing function

From the condition of sustained oscillation where is the linear element in one‐transcriptional module. From this we get one magnitude and one angle condition where and . This yields and

Calculation of the input bias

The simulation suggests that there is no unique input bias and amplitude pair (a, b) which gives rise to sustained oscillation. So, we need to fix the input bias. Considering the first transcription module in (11), the nonlinear part can be approximated to, . is the bias term in the approximation and the negative sign suggests phase for describing function approximation (Fig. 6). N is calculated from (15). Now, putting and , we get and Equating the constant terms we get which gives an equation of straight line, where every point corresponds to the condition of sustained oscillation. Now, from simulation also we can compute the value of when input bias a is varied which produces a family of curves with varying dependence of input bias a on input amplitude b. The intersection points correspond to the value of input bias. The simulation results are closest to the intersection points found by this method for a = b, when n = 2.

Describing function approximation with error estimate

The magnitude of the describing function approximation with the upper and lower error limits is shown in Fig. 6. We find the frequency of oscillation from (14) and this matches with simulation result. The describing function gain is calculated from (15) and corresponding input bias is calculated from (16). From the sustained oscillation condition, we get the nominal amplitude of oscillation. The upper limit of the error gives an upper limit of amplitude. The lower limit of error does not give a clear result as the error magnitude is larger than the nominal value. As we have taken same production and degradation constants for all three transcription module, we assume that the phase difference between each protein concentration profile is . The actual simulation of the repressilator system is compared with the oscillation profiles obtained from this analysis in Figs. 6 for different parameter values. We conclude that such as in the case of the Van der Pol oscillator, the use of the describing function technique with the simple error model adds to the limit cycle analysis by providing an estimate of the variation in the oscillation amplitude.

Discussion

Describing functions can be used to approximate a nonlinear input–output map with its linearisation. Here, we have adapted this method for investigations of biomolecular systems and presented the following three results. First, we used this technique to approximate representative input–output responses, both saturating and hysteretic, and mapped the dependence of the approximation on system parameters. Second, we estimated the approximation error, which was smaller than that of a standard linearisation, and theoretically developed a way to obtain an upper bound of this error. Third, we used the computed error estimate to augment the standard limit cycle prediction by providing a simple way to estimate range of oscillation amplitudes possible. These results should help to augment a framework for approximating biomolecular signalling systems with linearised versions.

It is interesting to note the additional insight that can be obtained by contrasting the describing function‐based linearisation with the standard linearisation. For example, in the high sensitivity regime of the enzymatic signalling system (Fig. 2), standard linearisation at the operating point would point to a dynamic response which is very slow and has a high gain. However, an analysis using a finite amplitude input, as in the describing function approximation here, adds insight to the fact that the slow response is present only if the amplitude is infinitesimal. For more realistic inputs, where the amplitude may be finite, the response is not significantly slower than in the low sensitivity regime (Fig. 2), though the gain is still higher. This provides a holistic understanding of this signalling system. In case of positive feedback systems, describing function approximation gives an account of the delay (Fig. 3) associated with it in the hysteretic regime. In the presence of multiple steady states, when standard linearisation is not well‐defined, describing function approximation gives a quantitative measure of this delay and can further aid in the design of biomolecular oscillators based on this.

An important task for the future is to obtain describing function approximations for different systems, both for building blocks such as analysed over here, as well as for analysis of larger systems obtained by combining such smaller systems. Furthermore, it may be useful to compare these representations with relatively recent experimental data on the response of such systems to periodic input [7, 8, 9, 10]. For example, we have computed describing function approximation for the input–output map of the galactose metabolism pathway in S. cerevisiae using a model described in [10] consisting of 18 differential equations. This is reported in supplementary material S2. We find that for this case also describing function serves as a better approximation than standard linearisation and can aid the analysis of frequency response data obtained experimentally. It may also be useful to analyse such complex networks by dividing them into multiple stages and replacing each with their linear approximations. Such a cascade network where each stage is a two component signalling system is analysed in supplementary material S3. We have shown how each stage can be approximated with a describing function‐based approximation and the overall approximation has lower error than standard linearisation. Similar approach can be adopted for larger biomolecular systems with the aim of enhancing ease of analysis. Finally, the use of random inputs to obtain such describing function approximations may help in treating biomolecular noise that can be potentially important [31].

A mathematical framework can help to understand the behaviour of large complex systems and as a tool for design. Additionally, an appreciation of the limitations of the framework undeniably aids these aims. Here, we have adapted the method of describing functions for approximating biomolecular systems and estimated the corresponding error. These can be also used, for example, in estimating the output response to other periodic inputs such as square waves (Fig. 7). To illustrate, we compare the output with a square wave to the one obtained by decomposing the square wave into three (fundamental, third and fifth harmonics) Fourier components and then adding the output response of these components as predicted by the describing function. The output achieved through these two approaches is reasonably similar. Finally, it should aid in developing a framework for analysis and design of larger, more complex biomolecular systems through systematic interconnections of smaller components.

Fig. 7

Application of describing function technique in estimating output to square wave input for the system in Fig. 1 for cycles/h. Square wave input is approximated up to the fifth harmonic and each harmonic is treated as a sinusoidal input for the use of describing function. Output to each harmonic summed together to get the approximate output and the results are compared

Supplementary Data 1

Click here for additional data file.