On the role of intrinsic noise on the response of the p53-Mdm2 module

Probably the simplest version, and the starting point of any development of the p53-Mdm2 module is the loop presented in figure 1. This is a feedback loop in which p53 acts as a transcription factor that regulates the expression of Mdm2 which, in turn, regulates the activity of p53. This is a well studied model [8, 9], whose mean field solution is characterized by damped oscillations. Sometimes, it is convenient to introduce in the model an intermediary species for the interaction between p53 and the Mdm2, usually a messenger mRNA that acts as the transcription factor of the Mdm2. We have shown before [19] that this intermediary is the responsible for the above mentioned dephasing between the oscillations of p53 and Mdm2 proteins. However, for simplicity we will assume here that it is in equilibrium, keeping the scheme of figure 1, as our basic block for the p53-Mdm2 module.

Zoom In Zoom Out Reset image size

We formalize the model starting from a stochastic approach considering the number of components of each species, p53 (${{n}_{\text{p}53}}$) and Mdm2 (${{n}_{\text{Mdm}2}}$) as the relevant variables of the problem.

From this point of view one faces the problem of solving the following Master equation:

$$\begin{eqnarray}&&\frac{\partial P(\vec{n},t)}{\partial t}= \sum\limits_{\vec{n}\neq\vec{n^{\prime}}} (P(\vec{n^{\prime}},t)W_{\vec{n^{\prime}}\rightarrow\vec{n}}-P(\vec{n},t)W_{\vec{n}\rightarrow\vec{n^{\prime}}})\end{eqnarray} \tag{ 1 }$$

where $P\left(\vec{n},t\right)$ is the probability that the system is in state $\vec{n}=\left({{n}_{\text{p}53}},{{n}_{\text{Mdm}2}}\right)$ at time t and ${{W}_{\vec{n}\to \overrightarrow{{{n}^{\prime}}}}}$ are transition rates from one state to another. The solution of equation (1) is an easy task only for very simple transition rates and one must usually use numerical methods, or resort to approximations. In this work we try both approaches and compare the results.

The simulations of the temporal evolution of the system can be obtained using the Gillespie algorithm [20]. As other Monte Carlo techniques this algorithm guarantees an exact solution, but it is time consuming and is not a good tool to explore the parameter space. Moreover, to understand the role of the parameters of the model it is important to have at least approximate solutions to the problem. Here, we use the Linear Noise Approximation (LNA) [21], that was first proposed in the context of chemical kinetics and have gain considerable attention in the last few years for modeling intracellular processes, [22–24], but also in more general contexts [25, 26]. In addition we compare the results of this approximation with experimental data from the literature. The mathematical details can be found in [21], and we summarize them in the appendix A for completeness. However, to simplify the reading, and to keep the consistency of the manuscript we present here its main results.

One of the difficulties to solve equation (1) is the discrete character of n_i. To deal with this, van Kampen [21] proposed the following approximation:

$$\begin{eqnarray}&&{{n}_{i}}=\Omega {{x}_{i}}+\sqrt{\Omega}{{\alpha}_{i}}\end{eqnarray} \tag{ 2 }$$

i.e. one approximates the number of molecules of species i by some concentration, plus some fluctuations proportional to the square root of the system size $\Omega$. Substituting (2) into (1) and developing the equation in order of $\frac{1}{\Omega}$ one obtains, at first order a set of differential equations for x_i, usually called the mean field or deterministic approximation:

$$\begin{eqnarray}&&\frac{\partial {{x}_{i}}}{\partial \tau}=\underset{j}{\mathop \sum}{{S}_{ij}}\widehat{{{W}_{j}}}\left(\vec{x}\Omega \right)\end{eqnarray} \tag{ 3 }$$

where S_ij is the stoichiometric matrix of the interactions (see the appendix A) for a detail derivation). At second order one can prove that the fluctuations ${{\alpha}_{i}}$ are Gaussian variables. However, to compare them with the experiments it is more useful to describe the evolution of these variables in terms of Langevin equations.

$$\begin{eqnarray}&&\frac{\partial {{\alpha}_{i}}}{\partial \tau}=\underset{k=1}{\overset{N}{\mathop \sum}}{{A}_{ik}}{{\alpha}_{k}}+{{\eta}_{i}}\left(\tau \right)\end{eqnarray} \tag{ 4 }$$

where ${{A}_{ik}}=\sum\nolimits_{j=1}^{M}{{S}_{ik}}\frac{\partial {{W}_{ij}}}{\partial {{x}_{k}}}$ and the ${{\eta}_{i}}$ has zero variance and correlation function $\langle {{\eta}_{i}}\left(\tau \right){{\eta}_{j}}\left({{\tau}^{\prime}}\right)\rangle ={{D}_{ik}}\delta \left(\tau -{{\tau}^{\prime}}\right)$ with ${{D}_{ik}}=\sum\nolimits_{j=1}^{M}{{S}_{ij}}{{S}_{kj}}{{W}_{j}}(\vec{x})$ (see the appendix A).

Equation (4) is a linear Langevin equation that can be easily solved in the Fourier representation

$$\begin{eqnarray}&&{{\hat{\alpha}}_{i}}(w)=\underset{k}{\mathop \sum}\Phi_{ik}^{-1}{{\hat{\eta}}_{k}}(w)\end{eqnarray} \tag{ 5 }$$

where ${{\Phi}_{ik}}=\text{i}w{{\delta}_{ik}}-{{A}_{ik}}$. From ${{\hat{\alpha}}_{i}}(w)$ one can get the power spectrum of the fluctuations as:

$$\begin{eqnarray}&&{{P}_{i}}(w)=\langle \vert {{\alpha}_{i}}(w){{\vert}^{2}}\rangle =\sum_j\sum_k \Phi_{ij}^-1(\omega)D_{jk}(\Phi^*)^{-1}_{ki}(\omega)\end{eqnarray} \tag{ 6 }$$

which, when studying fluctuations in single cells, is the quantity usually reported experimentally [9].

In the absence of any external signals p53, is activated only by random DNA damage. This is a common event during the cell life cycle. Spontaneous hydrolysis, collapsing replication forks and oxidative stress [27, 28], but also metabolic stress may produce p53 activation [5]. Within this context we study the basal response using the model in figure 1.

In stochastic terms the model is well described by the following set of expressions:

$$\begin{eqnarray}&&\text{p}53\overset{{{W}_{1}}}{\longrightarrow}\text{p}53+\text{Mdm}2\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\text{p}53\overset{{{W}_{2}}}{\longrightarrow}2\text{p}53\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\text{p}53\overset{{{W}_{3}}}{\longrightarrow}\Phi\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\text{Mdm}2\overset{{{W}_{4}}}{\longrightarrow}\Phi\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\text{p}53+\text{Mdm}2\overset{{{W}_{5}}}{\longrightarrow}\text{Mdm}2\end{eqnarray} \tag{ 11 }$$

where equations (7) and (8) reflect the activation of Mdm2 and the self-activation of p53. Equations (9) and (10) the decay of both species and equation (11) reflects the regulation of p53 by Mdm2. The species $\Phi$ reflects a vacuum state that is irrelevant in the analytical description of the problem, but is conveniently introduced to allow the free variation of the quantities of p53 and Mdm2 during the simulation keeping the system's size controlled [24]. The intuition is that if N is the total number of particles in the system, these are p53, Mdm2, and the number of molecules of the $\Phi$ species. In this sense, the dummy variable can be interpreted as the set of species in the cell that are connected with the process of interest, but not taken explicitly into account within the model.

The different W_i are transition probabilities proportional to the size of the system [21]. To fix concepts we work with the simplest possible transition probabilities, in this case:

$$\begin{eqnarray}&&{{W}_{1}}={{k}_{1}}\Omega {{x}_{\text{p}53}}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{W}_{2}}={{k}_{2}}\Omega {{x}_{\text{p}53}}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{W}_{3}}={{k}_{3}}\Omega {{x}_{\text{p}53}}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{{W}_{4}}={{k}_{4}}\Omega {{x}_{\text{Mdm}2}}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{{W}_{5}}={{k}_{5}}\Omega {{x}_{\text{p}53}}{{x}_{\text{Mdm}2}}\end{eqnarray} \tag{ 16 }$$

The transitions (W_{1, 2}), i.e. activation of Mdm2, and self-activation of p53 are proportional to the amount of p53 in the system. The degradation of Mdm2, (W₄) is proportional to the amount of Mdm2. The p53 is degraded following two mechanisms, a normal degradation (W₃) proportional to the quantity of p53 in the cell, and a Mdm2-mediated degradation (W₅), that is proportional to the quantities of p53 and of Mdm2. This is a clear simplification of the problem. However, one must first notice that non-linear approaches, although used with success in the modelling of gene regulatory networks are approximate expressions designed to describe a known biological phenomenology in systems governed by large populations, where the concept of 'concentration' is well defined [29, 30]. On the contrary we are dealing with problems where the number of molecules of the species involved is rather small, and therefore the concept itself of concentration is ill-defined. Of course, one could also write down transition rates that may reproduce in the proper limit of large number of molecules Hill-type functions, but it is not necessary true that this is a better assumption for small systems. Our linear hypothesis is nevertheless valid far from the saturation regime of the Hill function- that indeed is usually defined for very large concentrations. Therefore, following (3) the mean-field equations of this model read:

$$\begin{eqnarray}&&\frac{\text{d}{{x}_{\text{p}53}}}{\text{d}t}={{k}_{2}}{{x}_{\text{p}53}}-{{k}_{3}}{{x}_{\text{p}53}}-{{k}_{5}}{{x}_{\text{p}53}}{{x}_{\text{Mdm}2}}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\frac{\text{d}{{x}_{\text{Mdm}2}}}{\text{d}t}={{k}_{1}}{{x}_{\text{p}53}}-{{k}_{4}}{{x}_{\text{Mdm}2}}\end{eqnarray} \tag{ 18 }$$

The fixed points of these equations are: (0, 0) and $\left(\frac{{{k}_{4}}{{k}_{2}}}{{{k}_{1}}{{k}_{5}}},\frac{{{k}_{2}}}{{{k}_{5}}}\right)$ where to simplify notation we considered $k_{2}^{\prime}={{k}_{2}}-{{k}_{3}}$ and relabeled it as k₂. Then, the non trivial solution has biological sense only if the rate of activation of p53 is larger than its degradation rate (k₂  >  0).

The stability of these fixed points is defined by the following Lyapunov exponents:

$$\begin{eqnarray}&&(0,0)\to \vec{\lambda}=\left\{\begin{array}{r} {{k}_{2}}\\ -{{k}_{4}} \end{array}\right\} \end{eqnarray} \tag{ 19 }$$

and:

$$\begin{eqnarray}&&\left(x_{\text{p}53}^{*},x_{\text{Mdm}2}^{*}\right)\to \vec{\lambda}=\left\{\begin{array}{c} \frac{1}{2}\left(-{{k}_{4}}+\sqrt{k_{4}^{2}-4{{k}_{2}}{{k}_{4}}}\right)\\ \frac{1}{2}\left(-{{k}_{4}}-\sqrt{k_{4}^{2}-4{{k}_{2}}{{k}_{4}}}\right) \end{array}\right\}\end{eqnarray} \tag{ 20 }$$

These results guarantee that if the non-trivial solution exists (k₄  > 0 and k₂  >  0), then it is stable. Moreover, if ${{k}_{4}}<4{{k}_{2}}$, i.e. if the activation of p53 is large enough, it displays damped oscillations. This is the region of parameters relevant from the biological point of view.

In figure 2 we show the comparison between the numerical solution of these equations with the results of the Gillespie algorithm for the same values of the parameters. Note that, while the deterministic solution is over-damped, the stochastic solution displays persistent oscillations around the fixed point.

Zoom In Zoom Out Reset image size

To characterize these oscillations and in order to compare our model with the experimental results reported in [9], we computed the power spectrum of the results of the Gillespie simulations over 200 realizations and made their geometrical average. The resulting power spectrum of the oscillations of p53 extracted from the simulations is presented (symbols) in figure 3. In the same figure we also show (line) the analytical results predicted by the LNA (see equation (6) and the appendix A). The inset presents the tail of the power spectrum in a logarithmic scale. The linear dependency of this tail with the inverse of the frequency is also a prediction of the LNA.

Zoom In Zoom Out Reset image size

The existence of a maximum in the power spectrum of the fluctuations in figure 3 is associated to the presence of a characteristic frequency for these fluctuations. In an infinite system, these oscillations die out, but the stochasticity due to the finite size of the system acts as a white noise, see equation (4), that interacts with the proper frequency of the system producing a resonance-like effect. This characteristic frequency is connected with the presence of an imaginary part in the eigenvalues of the system of differential equations and is a phenomena that has been very well characterized before in other models, see for example, [22–26]. In our specific case it is easy to show (see the appendix A), that the frequency of the oscillations only depends on the effective activation rate of p53 (${{k}_{2}}-{{k}_{3}}$) and the decay rate of Mdm2(k₄). On the contrary, the activation rate of Mdm2 by p53 (k₁) and the regulation of the p53 by Mdm2(k₅) although unavoidable to build the feed-back loop in the model are not connected with the characteristic frequency of the oscillations, suggesting a possible robust mechanism of interaction with other species.

The comparison with the experimental data is simplified with the help of the previous theoretical predictions that set a relation between some parameters of the model and the frequency of oscillations. This in turn can be extracted from the experimental data [9]. In this case for example we know (see the appendix A) that the frequency of oscillations is defined through the expression:

$$\begin{eqnarray}&&\omega =2\pi f\approx \sqrt{\left({{k}_{2}}-{{k}_{3}}\right){{k}_{4}}}\end{eqnarray} \tag{ 21 }$$

and from the experimental data [9] that f   =  0.08 h⁻¹, it is simply then to choose values of ${{k}_{2}},{{k}_{3}}$ and k₄ that satisfy this condition and are consistent with previous predictions in the literature.

In fact, figure 3 compares also very well with figure 6 in [9] where the authors studied the basal response of the system and presented experimental results for the Fourier transform of the fluctuations. In that work the authors model their results adding some noise ad-hoc to mean field equations that are similar to ours. On the contrary, in our approach, the oscillations arise naturally from the finite size of the populations involved and reproduce very well the value reported in the experiments f   =  0.08 h⁻¹.

The rest of the parameters chosen have little influence in the main characteristic of the curve and as can be seen in the table 1 they compare very well with previous reports in the literature. Only k₅, the parameter reflecting the Mdm2-depending degradation rate of p53 is clearly underestimated in comparison with the value reported in [8]. But one must notice that, on one hand, their activation rate of p53 (k₂) is larger than ours, and on the other, that our model includes and Mdm2-independent degradation rate of p53 (k₃) that is absent in the model V of reference [8]. Therefore this underestimation of k₅ is not surprising at all.

Table 1. Comparison between our predictions for the values of the parameters of the model and similar estimates in the literature (Model V in [8]).

Parameter h−1	This work	Literature
k1	.99	1.5 $\pm$60% [8]
k2	1.0	2.0 $\pm$25% [8]
k3	0.44	—,— [8]
k4	.49	0.9 $\pm$30% [8]
k5	1.05	3.7 $\pm$50% [8]

We already mentioned that p53 is one of the central proteins in the cellular response to stress. One of the most typical stresses is DNA double strand break (DSB), which could be generated, for example, by gamma irradiation. The DSBs produces the activation of the upstream kinase ATM, that in turn induces the phosphorylation (activation) of p53. The p53, as already mentioned, then works as a transcription factor of several target genes. These genes activate different responses, ranging from DNA repair to apoptosis [10].

It is then natural, when modelling this kind of stress, to incorporate ATM as a new dynamical species. In figure 4 we present a sketch of this model. From the technical point of view, we must add to the dynamics governing the basal state the following new transitions:

$$\begin{eqnarray}&&\text{ATM}\overset{{{W}_{6}}}{\longrightarrow}\text{ATM}+p53\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\text{ATM}+p53\overset{{{W}_{7}}}{\longrightarrow}p53\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\text{ATM}\overset{{{W}_{8}}}{\longrightarrow}\Phi\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\text{ATM}\overset{{{W}_{9}}}{\longrightarrow}2\text{ATM}\end{eqnarray} \tag{ 25 }$$

Zoom In Zoom Out Reset image size

These transitions reflect the activation of p53 by the ATM, (W₆), the regulation of the ATM by p53, (W₇), the natural degradation of the ATM, (W₈), and a self-reinforcement field for the ATM, (W₉) that mimics the persistence of the damage. It is important however to take in mind that W₇ summarized a complex set of processes that after activation by p53 let to the DNA repair. As in the previous section we assumed that they take the following simple form:

$$\begin{eqnarray}&&{{W}_{6}}={{k}_{6}}\Omega {{x}_{\text{ATM}}}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{{W}_{7}}={{k}_{7}}\Omega {{x}_{\text{p}53}}{{x}_{\text{ATM}}}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{{W}_{8}}={{k}_{8}}\Omega {{x}_{\text{ATM}}}\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{{W}_{9}}={{k}_{9}}\Omega {{x}_{\text{ATM}}}\end{eqnarray} \tag{ 29 }$$

With these definitions of the transition rules, we may write the mean field equations of the model:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}\frac{\text{d}{{x}_{\text{p}53}}}{\text{d}t}={{k}_{2}}{{x}_{\text{p}53}}-{{k}_{3}}{{x}_{\text{p}53}}\\ -{{k}_{5}}{{x}_{\text{p}53}}{{x}_{\text{Mdm}2}}+{{k}_{6}}{{x}_{\text{ATM}}}\\ \frac{\text{d}{{x}_{\text{Mdm}2}}}{\text{d}t}={{k}_{1}}{{x}_{\text{p}53}}-{{k}_{4}}{{x}_{\text{Mdm}2}}\\ \frac{\text{d}{{x}_{\text{ATM}}}}{\text{d}t}={{k}_{9}}{{x}_{\text{ATM}}}-{{k}_{8}}{{x}_{\text{ATM}}}-{{k}_{7}}{{x}_{\text{ATM}}}{{x}_{\text{p}53}}\end{array}\end{eqnarray} \tag{ 30 }$$

These equations are similar to the ones in (18), but now the role of ATM on the activation of p53 results in the addition of a new term to the first equation. Moreover, ATM is also a new dynamical variable whose evolution must be taken into account.

This system of equations has three different fixed points. They are, (0, 0, 0), $\left(\frac{{{k}_{2}}{{k}_{4}}}{{{k}_{1}}{{k}_{5}}},\frac{{{k}_{2}}}{{{k}_{4}}},0\right)$ and $\left(\frac{{{k}_{9}}}{{{k}_{7}}},\frac{{{k}_{9}}{{k}_{1}}}{{{k}_{7}}{{k}_{4}}},\frac{{{k}_{9}}}{{{k}_{7}}{{k}_{6}}}\left(\frac{{{k}_{9}}{{k}_{1}}{{k}_{5}}}{{{k}_{7}}{{k}_{4}}}-{{k}_{2}}\right)\right)$ where for simplicity we re-labelled again $k_{2}^{\prime}={{k}_{2}}-{{k}_{3}}$, and $k_{9}^{\prime}={{k}_{9}}-{{k}_{8}}$, as k₂ and k₉ respectively.

The Jacobian matrix necessary to study their stability takes the general form:

$$\begin{eqnarray}J=\left[\begin{array}{c} {{k}_{2}}-{{k}_{5}}x_{\text{Mdm}2}^{*} & -{{k}_{5}}x_{\text{p}53}^{*} & {{k}_{7}}\\ {{k}_{1}} & -{{k}_{4}} & 0\\ -{{k}_{7}}x_{\text{ATM}}^{*} & 0 & 0 \end{array}\right]\end{eqnarray}$$

that evaluated in the fixed point of interest $\left(\frac{{{k}_{9}}}{{{k}_{7}}},\frac{{{k}_{9}}{{k}_{1}}}{{{k}_{7}}{{k}_{4}}},\frac{{{k}_{9}}}{{{k}_{7}}{{k}_{6}}}\left(\frac{{{k}_{9}}{{k}_{1}}{{k}_{5}}}{{{k}_{7}}{{k}_{4}}}-{{k}_{2}}\right)\right)$ transforms into the following secular equation:

$$\begin{eqnarray}&&{{\lambda}^{3}}+{{\lambda}^{2}}\left(c-{{k}_{4}}\right)+\lambda \left(\frac{{{k}_{9}}{{k}_{1}}{{k}_{5}}}{{{k}_{7}}}-{{k}_{4}}c-{{k}_{9}}c\right)-{{k}_{9}}{{k}_{4}}c=0\end{eqnarray} \tag{ 31 }$$

where $c={{k}_{2}}-\frac{{{k}_{9}}{{k}_{1}}{{k}_{5}}}{{{k}_{7}}{{k}_{4}}}$.

This secular equation must be studied numerically and its solution defines the stability of the fixed points. We have found that depending on the parameters, three different regions are well defined. In Region I the solution $\{x_{\text{p}53}^{*},x_{\text{Mdm}2}^{*},x_{\text{ATM}}^{*}\}$ may exist, is stable and the system reaches the fixed point with damped oscillations. In this region also the solution $\{x_{\text{p}53}^{*},x_{\text{Mdm}2}^{*},0\}$ is stable if k₂  >  0. Depending on the initial condition, the system is attracted to one or another fixed point. In Region II the solution $\{x_{\text{p}53}^{*},x_{\text{Mdm}2}^{*},x_{\text{ATM}}^{*}\}$ has biological sense and is stable, but the fixed point is reached without damped oscillations. In this region the solution $\{x_{\text{p}53}^{*},x_{\text{Mdm}2}^{*},0\}$ does not exist. Finally, in Region III only the solution $\{x_{\text{p}53}^{*},x_{\text{Mdm}2}^{*},0\}$ is stable.

Based on the results of the previous section we concentrate our attention in Region I, where the mean field fixed point $\{x_{\text{p}53}^{*},x_{\text{Mdm}2}^{*},x_{\text{ATM}}^{*}\}$ is reached quickly through over-damped oscillations. In figure 5 we show the numerical solution of equation (29) and the results of the Gillespie algorithm using the parameters $\vec{k}=\left(0.99,1,0.44,0.69,0.85,0.5,0.5,0.1,0.4\right){{h}^{-1}}$. In the left panel of the figure we plot the first 30 h of the simulations, and compare them with the mean-field solutions. It is clear that while the over-damped oscillations die out very fast, the stochastic fluctuations are sustained by a large period of time. This is represented in the right panel in the same figure where these fluctuations are shown for up to 110 h.

Zoom In Zoom Out Reset image size

To characterize these oscillations we use again the power spectrum of the Gillespie results and compare them with the predictions from the LNA. In figure 6 we show with points this power spectrum and with a continuous line the predictions from the LNA. As in the previous section, in the inset we plot the tail power spectrum in a logarithmic scale, showing again a linear dependency with $1/\omega$. Notice the good coincidence between the numerical simulations and the LNA.

Zoom In Zoom Out Reset image size

We may also compare our results with figure 3 in [9]. Again, our model may reproduce qualitatively and quantitatively the experimental observations: the existence of a characteristic frequency (f   =  0.14 h⁻¹), compatible with the experimental results and larger that in the basal state, and a divergence close to zero frequency. The parameters were intentionally chosen to resemble the experiments. As before, the characteristic frequency of the oscillations is essentially defined by ${{k}_{2}}-{{k}_{3}}$ and k₄. We fit the expected experimental frequency (f   =  .14 h⁻¹) [9], increasing (decreasing) k₄ (k₅) with respect to the model in the basal state to resemble the experimentally verified fact that under DNA damage the degradation rate of Mdm2 increases and the Mdm2-mediated degradation rate of p53 decreases [31]. Moreover, to reproduce the divergence at low frequency of the model we noticed that it is associated with the parameters characterizing the new feed-back loop. It can be interpreted looking to the low frequency analytical expression of the power spectrum derived from the LNA. It has the non-trivial form:

$$\begin{eqnarray}&&{{P}_{\text{p}53}}\left(w\to 0\right)\sim \frac{2{{k}_{5}}{{k}_{7}}}{{{k}_{1}}{{k}_{6}}{{k}_{8}}-{{k}_{2}}{{k}_{5}}{{k}_{7}}}\end{eqnarray} \tag{ 32 }$$

that nonetheless guarantees, if the denominator is small enough, a large value for the power spectrum. With this expression, the numerical results of the previous section and previous estimates from the literature, it was easy to choose values for the parameters that reproduce the divergence of the experimental data and the correct frequency of resonance (see table 2). However, we want to stress that there is still a lot of freedom in the parameter selection and our interest here was not to find the best ones, but just to show that this approach provides a straightforward path to obtain a quantitative comparison between our theoretical model and the experimental data. It must be noticed that excluding the parameters that were fixed to reproduce the divergence at low frequency and the peak in the power spectrum, variations within a 10% of the values reported in table 2 reproduce curves qualitatively similar to figure 6. All these values compare very well with previous theoretical estimates of equivalent parameters in similar models, as presented also in table 2.

Table 2. Comparison between our predictions for the values of the parameters of the model and similar estimates in the literature, Model VI in [8] and a cross-regulated linear model in [9].

Parameter h−1	This work	Literature
k1	0.99	$0.29\pm .03$% [9], 1.20 [8]
k2	1	$2\pm .5$% [9], 0.8 [8]
k3	.44	$.45\pm .03$% [9], $1.40\pm 20$% [8]
k4	.69	$.28\pm .02$% [9]
k5	.85	$.55\pm .05$% [9]
k6	.5	${{k}_{6}}{{k}_{7}}=.65\pm .05$% [9]
k7	.5	${{k}_{6}}{{k}_{7}}=.65\pm .05$% [9]
k8	.4

Another important pathway for the activation of p53, comes from the activation of the p14^ARF protein. In the presence of oncogenes this protein is activated and joins Mdm2 accelerating its degradation [32]. This process leads to the increase of p53 level in the cell. Unfortunately, as far as we know there is no experimental evidence for the single cell dynamics of p53 and Mdm2 in the presence of oncogenes, therefore some of the parameters used in this section lack of experimental bias. Nevertheless we will show that independently of the values used, the power spectrum of the fluctuations is qualitatively different from the one obtained when p53 interacts directly with ATM.

The model of oncogenic response is sketched in figure 7. To the standard interactions present in the basal model we add the regulation of p14^ARF to Mdm2 and a regulation of p14^ARF by p53. In terms of stochastic transitions we have now:

$$\begin{eqnarray}&&\text{Mdm}2+p{{14}^{\text{ARF}}}\overset{{{W}_{6}}}{\longrightarrow}p{{14}^{\text{ARF}}}\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&p{{14}^{\text{ARF}}}+p53\overset{{{W}_{7}}}{\longrightarrow}p53\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&p{{14}^{\text{ARF}}}\overset{{{W}_{8}}}{\longrightarrow}\Phi\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&p{{14}^{\text{ARF}}}\overset{{{W}_{9}}}{\longrightarrow}2p{{14}^{\text{ARF}}}\end{eqnarray} \tag{ 36 }$$

and defining the transition rates as:

$$\begin{eqnarray}&&{{W}_{6}}={{k}_{6}}\Omega {{x}_{\text{p}14}}{{x}_{\text{Mdm}2}}\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{{W}_{7}}={{k}_{7}}\Omega {{x}_{\text{p}53}}{{x}_{\text{p}14}}\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{{W}_{8}}={{k}_{8}}\Omega {{x}_{\text{p}14}}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{{W}_{9}}={{k}_{9}}\Omega {{x}_{\text{p}14}}\end{eqnarray} \tag{ 40 }$$

we obtain the following system of equations for the dynamic of the model:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}\frac{\text{d}{{x}_{\text{p}53}}}{\text{d}t}={{k}_{2}}{{x}_{\text{p}53}}-{{k}_{3}}{{x}_{\text{p}53}}-{{k}_{5}}{{x}_{\text{p}53}}{{x}_{\text{Mdm}2}}\\ \frac{\text{d}{{x}_{\text{Mdm}2}}}{\text{d}t}={{k}_{1}}{{x}_{\text{p}53}}-{{k}_{4}}{{x}_{\text{Mdm}2}}-{{k}_{6}}{{x}_{\text{p}14}}{{x}_{\text{Mdm}2}}\\ \frac{\text{d}{{x}_{\text{p}14}}}{\text{d}t}={{k}_{9}}{{x}_{\text{p}14}}-{{k}_{8}}{{x}_{\text{p}14}}-{{k}_{7}}{{x}_{\text{p}14}}{{x}_{\text{p}53}}\end{array}\end{eqnarray} \tag{ 41 }$$

which again depending on the parameters may show three fixed points: (0, 0, 0), $\left(\frac{{{k}_{2}}{{k}_{4}}}{{{k}_{1}}{{k}_{5}}},\frac{{{k}_{1}}}{{{k}_{5}}},0\right)$ and $\left(\frac{{{k}_{9}}}{{{k}_{7}}},\frac{{{k}_{2}}}{{{k}_{5}}},\frac{1}{{{k}_{2}}{{k}_{6}}}\left(\frac{{{k}_{9}}{{k}_{1}}{{k}_{5}}}{{{k}_{7}}}-{{k}_{2}}{{k}_{4}}\right)\right)$. To simplify the notation we have used ${{k}_{2}}={{k}_{2}}-{{k}_{3}}$ and ${{k}_{9}}={{k}_{9}}-{{k}_{7}}$.

Zoom In Zoom Out Reset image size

The Jacobian matrix near these fixed points takes the form:

$$\begin{eqnarray}J=\left[\begin{array}{c} {{k}_{2}}-{{k}_{5}}x_{\text{Mdm}2}^{*} & -{{k}_{5}}x_{\text{p}53}^{*} & 0\\ {{k}_{1}} & -{{k}_{4}} & -{{k}_{6}}x_{\text{Mdm}2}^{*}\\ -{{k}_{2}}x_{\text{p}14}^{*} & 0 & {{k}_{9}} \end{array}\right]\end{eqnarray} \tag{ 42 }$$

that evaluated in the fixed point $\{x_{\text{p}53}^{*},x_{\text{Mdm}2}^{*},x_{\text{p}14}^{*}\}$ leads to the following secular equation:

$$\begin{eqnarray}&&{{\lambda}^{3}}+{{\lambda}^{2}}{{k}_{4}}+\lambda \left(\frac{{{k}_{9}}{{k}_{1}}{{k}_{5}}}{{{k}_{7}}}\right)+\frac{k_{9}^{2}{{k}_{5}}{{k}_{1}}}{{{k}_{7}}}-{{k}_{2}}{{k}_{9}}{{k}_{4}}=0\end{eqnarray} \tag{ 43 }$$

which defines a phase diagram qualitatively similar to the one discussed in the preceding section. Three regions, but only one of them, with solution $\{x_{\text{p}53}^{*},x_{\text{Mdm}2}^{*},x_{\text{p}14}^{*}\}\ne \vec{0}$ displays damped oscillations. Our previous results suggest that this is the relevant region for our purposes. In this region we used $\vec{k}=\left(1.0,1.0,0.44,0.69,0.85,0.5,0.5,0.4\right)$ as parameters to obtain damped oscillations in the mean field solution of the problem and sustained oscillations through the Gillespie algorithm. The results are shown in figure 8.

Zoom In Zoom Out Reset image size

In figure 9 we show the power spectrum obtained using the LNA and the one obtained from the simulations. Also for this model, the LNA approximation correctly describes the results of the simulation. As before a clear peak representative of the feed-back loop between p53 and Mdm2 is observed and also the exponential decay of the tail of the power spectrum with $1/\omega$. However, there is no evidence of a divergence at small frequencies.

Zoom In Zoom Out Reset image size

This absence of the divergence is independent of the parameters chosen. It can be understood looking again to the structure of the power spectrum predicted by the LNA at small frequencies. For this model:

$$\begin{eqnarray}&&{{P}_{\text{p}53}}\left(w\to 0\right)=\frac{2{{k}_{5}}\left({{k}_{1}}+{{k}_{2}}\right)}{k_{1}^{2}{{k}_{6}}}\end{eqnarray} \tag{ 44 }$$

that indicates that only when k₁ or k₆ are zero, i.e. only when the regulation of Mdm2 by p53 or the interaction between Mdm2 and p14^ARF disappears we have a diverging power spectrum at small frequencies. But in this case our model of interactions breaks down.

The quantities characterizing the number of particles in our problem are discrete. The deviation of this number from the expected mean value produces fluctuations usually called intrinsic noise because it is not caused by any external source. When the system is small enough, this intrinsic noise is not negligible, and a mean field approximation can not be taken.

In these cases, the correct way to describe the system is to solve the Master Equation that governs it:

$$\begin{eqnarray}&&\frac{\partial P(\vec{n},t)}{\partial t}= \sum\limits_{\vec{n}\neq\vec{n^{\prime}}} (P(\vec{n^{\prime}},t)W_{\vec{n^{\prime}}\rightarrow\vec{n}}-P(\vec{n},t)W_{\vec{n}\rightarrow\vec{n^{\prime}}}) .\end{eqnarray} \tag{ A.1 }$$

This is an equation for the change in time of the probability $P\left(\vec{n},t\right)$ of being in the state $\vec{n}$ at the time t, and it rests upon the Markovian hypothesis. ${{W}_{\vec{n}\to \overrightarrow{{{n}^{\prime}}}}}$ are transition rates from states $\vec{n}$ to $\overrightarrow{{{n}^{\prime}}}$.

An important difficulty in solving the Master equation arises from the discreteness of the variables involved ($\vec{n}$). Its solution is obtained by analytical methods in very lucky situations. An approximation scheme that allows to get some information from it, is the linear noise approximation (LNA).

This approximation assumes the possibility of studying the number n_i of integrands of the ith population as a main part $\Omega {{x}_{i}}$ that represents the mean value of the population size, plus a small deviation $\sqrt{\Omega}{{\alpha}_{i}}$ that represents the fluctuations around the mean value, and that is proportional to the square root of the system size $\Omega$. In this way:

$$\begin{eqnarray}&&{{n}_{i}}=\Omega {{x}_{i}}+\sqrt{\Omega}{{\alpha}_{i}}.\end{eqnarray} \tag{ A.2 }$$

As every process that takes place in the system causes variations of one or more integrands in the populations of the species involved, it is useful to write the master equation in terms of these discrete fluctuations. For this purpose let us define the step operator as:

$$\begin{eqnarray}&&E_{i}^{k}f\left(\ldots ,{{n}_{i}},\ldots \right)=f\left(\ldots ,{{n}_{i}}+k,\ldots \right).\end{eqnarray} \tag{ A.3 }$$

In terms of this operator the master equation reads:

$$\begin{eqnarray}&&\frac{\text{d}P\left(\vec{n},t\right)}{\text{d}t}=\underset{j=1}{\overset{M}{\mathop \sum}}\left[\left(\underset{i=1}{\overset{N}{\mathop \prod}}E_{i}^{-{{S}_{ij}}}\right)-1\right]{{W}_{j}}\left(\vec{n},\Omega \right)P\left(\vec{n},t\right),\end{eqnarray} \tag{ A.4 }$$

where W_j is the microscopic reaction rate for the jth process and S_ij is the element (i, j) of the stoichiometry matrix $\mathbf{S}$, the magnitude of the change in the ith population when the jth process takes place. The total number of different populations that form the system is N, and M is the total number of different kinds of processes that can take place in the system.

In terms of the fluctuations, which are continuous variables, the step operator takes the following form, more suitable for an analytic treatment

$$\begin{eqnarray}&&{{E}_{i}}^{k}\approx 1+\frac{k}{\sqrt{\Omega}}\frac{\partial}{\partial {{\alpha}_{i}}}+\frac{{{k}^{2}}}{2\Omega}\frac{{{\partial}^{2}}}{\partial {{\alpha}_{i}}^{2}}+\ldots \end{eqnarray} \tag{ A.5 }$$

in this way, we can work with a continuous operator, and still not lose the stochastic features of these processes. With it

$$\begin{eqnarray}&&\begin{array}{*{20}{r}}\underset{i=1}{\overset{N}{\mathop \prod}}E_{i}^{-{{S}_{ij}}}{{f}_{j}}\left(\vec{\alpha},\Omega \right)=\left[1-{{\Omega}^{-\frac{1}{2}}}\underset{i=1}{\overset{N}{\mathop \sum}}{{S}_{ij}}\frac{\partial}{\partial {{\alpha}_{i}}}\right]{{f}_{j}}\left(\vec{\alpha},\Omega \right)\\ +\frac{1}{2}\left[{{\Omega}^{-1}}\underset{i=1}{\overset{N}{\mathop \sum}}\underset{k=1}{\overset{N}{\mathop \sum}}{{S}_{ij}}{{S}_{kj}}\frac{\partial}{\partial {{\alpha}_{i}}}\frac{\partial}{\partial {{\alpha}_{k}}}+\ldots \right]{{f}_{j}}\left(\vec{\alpha},\Omega \right)\end{array}\end{eqnarray} \tag{ A.6 }$$

Now we are going to replace the probability $P\left(\vec{n},t\right)$ by the probability $\Pi(\vec{\alpha},t)$ of being $\sqrt{\Omega}\vec{\alpha}$ away from the mean field prediction

$$\begin{eqnarray}&&P\left(\vec{n},t\right)=P\left(\Omega \vec{x}+\sqrt{\Omega}\vec{\alpha}\right)=\Pi(\vec{\alpha},t)\end{eqnarray} \tag{ A.7 }$$

We know that:

$$\begin{eqnarray}&&\begin{array}{*{20}{r}}\frac{\partial \Pi(\vec{\alpha},t)}{\partial t}=\frac{\partial P\left(\vec{n},t\right)}{\partial t}+\Omega \underset{i=1}{\overset{N}{\mathop \sum}}\frac{\partial P\left(\vec{n},t\right)}{\partial {{n}_{i}}}\frac{\partial {{x}_{i}}}{\partial t}\;\;\;\;\;\;\text{and}\\ \;\;\;\;\;\;\frac{\partial P\left(\vec{n},t\right)}{\partial {{n}_{i}}}=\frac{1}{\sqrt{\Omega}}\frac{\partial \Pi(\vec{\alpha},t)}{\partial {{\alpha}_{i}}},\end{array}\end{eqnarray}$$

with this results we can obtain the relation

$$\begin{eqnarray}&&\frac{\partial \Pi(\vec{\alpha},t)}{\partial t}-{{\Omega}^{\frac{1}{2}}}\underset{i=1}{\overset{N}{\mathop \sum}}\frac{\partial \Pi(\vec{\alpha},t)}{\partial {{\alpha}_{i}}}\frac{\partial {{x}_{i}}}{\partial t}=\frac{\partial P\left(\vec{n},t\right)}{\partial t}.\end{eqnarray} \tag{ A.8 }$$

If now we expand the transition probabilities per unit time around $\vec{x}$

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{W}_{j}}(\vec{\alpha},\Omega )={{W}_{j}}\left(\vec{x}+{{\Omega}^{-\frac{1}{2}}}\vec{\alpha}\right)\simeq {{W}_{j}}(\vec{x})\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+{{\Omega}^{-\frac{1}{2}}}\underset{i=1}{\overset{N}{\mathop \sum}}\frac{\partial {{W}_{j}}(\vec{x})}{\partial {{x}_{i}}}{{\alpha}_{i}}+o({{\Omega}^{-1}})\end{array}\end{eqnarray} \tag{ A.9 }$$

Substituting (A.6) and (A.9) in the master equation (A.4) and using (A.7), after grouping the coefficients of the similar powers of $\Omega$ we find:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}\frac{\partial P\left(\vec{n},t\right)}{\partial t}={{\Omega}^{-1}}\left(\underset{j=1}{\overset{M}{\mathop \sum}}{{W}_{j}}(\vec{x})\underset{i,k=1}{\overset{N}{\mathop \sum}}{{S}_{ij}}{{S}_{kj}}\frac{{{\partial}^{2}}\Pi(\vec{\alpha},t)}{\partial {{\alpha}_{i}}\partial {{\alpha}_{k}}}\right)\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-{{\Omega}^{-1}}\left(\underset{j=1}{\overset{M}{\mathop \sum}}\underset{i,k=1}{\overset{N}{\mathop \sum}}{{S}_{ij}}\frac{\partial \left[\Pi(\vec{\alpha},t){{\alpha}_{k}}\right]}{\partial {{\alpha}_{i}}}\frac{\partial {{W}_{j}}(\vec{x})}{\partial {{x}_{k}}}\right)\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+{{\Omega}^{-\frac{1}{2}}}\left(-\underset{j=1}{\overset{M}{\mathop \sum}}\underset{i=1}{\overset{N}{\mathop \sum}}{{S}_{ij}}{{W}_{j}}(\vec{x})\frac{\partial \Pi(\vec{\alpha},t)}{\partial {{\alpha}_{i}}}\right)+o\left({{\Omega}^{-\frac{3}{2}}}\right)\end{array}\end{eqnarray}$$

If now we compare the coefficients of the similar powers of $\Omega$ in equations (A.10) and (A.8), where the unit of time has been changed from t to $\tau =\frac{t}{\Omega}$ we find for the leading order:

$$\begin{eqnarray}&&\frac{\partial {{x}_{i}}}{\partial \tau}=\underset{j=1}{\overset{M}{\mathop \sum}}{{S}_{ij}}{{W}_{j}}(\vec{x}).\end{eqnarray} \tag{ A.10 }$$

This corresponds to the deterministic rate equations that are often used to describe the dynamics of such systems at a mean-field level. It has been discussed already in the subsection referring to mean field equations and stability analysis.

For the coefficients of ${{\Omega}^{-1}}$ the relation reads:

$$\begin{eqnarray}&&\begin{array}{*{20}{r}}\frac{\partial \Pi\left(\alpha ,\tau \right)}{\partial \tau}=-\underset{i,k=1}{\overset{N}{\mathop \sum}}{{A}_{ik}}\frac{\partial \left[\Pi\left(\alpha ,t\right){{\alpha}_{k}}\right]}{\partial {{\alpha}_{i}}}\\ +\frac{1}{2}\underset{i,k=1}{\overset{N}{\mathop \sum}}{{D}_{ik}}\frac{{{\partial}^{2}}\Pi\left(\alpha ,t\right)}{\partial {{\alpha}_{i}}\partial {{\alpha}_{k}}},\end{array}\end{eqnarray} \tag{ A.11 }$$

where

$$\begin{eqnarray}&&{{A}_{ik}}=\underset{j=1}{\overset{M}{\mathop \sum}}{{S}_{ij}}\frac{\partial {{W}_{j}}(\vec{x})}{\partial {{x}_{k}}}\end{eqnarray} \tag{ A.12 }$$

is the jacobian matrix of the system, and

$$\begin{eqnarray}&&{{D}_{ik}}=\underset{j=1}{\overset{M}{\mathop \sum}}{{S}_{ij}}{{S}_{kj}}{{W}_{j}}(\vec{x}).\end{eqnarray} \tag{ A.13 }$$

Equation (A.11) is a Fokker-Plank equation for the probability for the system to have the deviation $\sqrt{\Omega}\vec{\alpha}$ around the mean-field prediction. This equation can be written in a completely equivalent formulation more benevolent to investigation using Fourier transforms. The problem can then be formulated as the set of stochastic differential equations of the Langevin type:

$$\begin{eqnarray}&&\frac{\text{d}{{\alpha}_{i}}}{\text{d}\tau}=\underset{k=1}{\overset{N}{\mathop \sum}}{{A}_{ik}}{{\alpha}_{k}}+{{\eta}_{i}}(t)\end{eqnarray} \tag{ A.14 }$$

where ${{\eta}_{i}}\left(\tau \right)$ is a Gaussian noise with 0 mean and correlation given by $\langle {{\eta}_{i}}\left(\tau \right){{\eta}_{k}}\left({{\tau}^{\prime}}\right)\rangle ={{D}_{ik}}\delta \left(\tau -{{\tau}^{\prime}}\right)$.

This a linear Langevin equation that can be easily solved in the Fourier representation

$$\begin{eqnarray}&&{{\hat{\alpha}}_{i}}(w)=\underset{k}{\mathop \sum}\Phi_{ik}^{-1}{{\hat{\eta}}_{k}}(w)\end{eqnarray} \tag{ A.15 }$$

where ${{\Phi}_{ik}}=iw{{\delta}_{ik}}-{{A}_{ik}}$. From ${{\hat{\alpha}}_{i}}(w)$ one can get the the power spectrum of the fluctuations as:

$$\begin{eqnarray}&&{{P}_{i}}(w)=\langle \vert {{\tilde{\alpha}}_{i}}(w){{\vert}^{2}}\rangle =\underset{j}{\mathop \sum}\underset{k}{\mathop \sum}\Phi_{ij}^{-1*}\Phi_{ik}^{-1}{{D}_{jk}}\end{eqnarray} \tag{ A.16 }$$

Now we are going to analyze the general expression for the power spectrum (A.16) for the basal state. The analytics for the other two models develops following the same steps, with a more cumbersome Algebra

$$\begin{eqnarray}&&{{P}_{i}}(\omega)\propto \underset{j=1}{\overset{N}{\mathop \sum}}\underset{k=1}{\overset{N}{\mathop \sum}}\Phi_{ij}^{-1}(\omega){{D}_{jk}}\left({{\Phi}^{\dagger}}\right)_{ki}^{-1}(\omega)\end{eqnarray} \tag{ A.17 }$$

Since $\Phi=\text{i}\omega {{\delta}_{ik}}-{{A}_{ik}}$ and since A and D are independent of ω, ${{P}_{i}}(\omega)$ is a fraction, the numerator and the denominator have both a polynomial structure of order 2N. The explicit form of the denominator is: $\vert \text{det} \Phi(\omega){{\vert}^{2}}$.

Then, we expect that the enhancement of the fluctuations would be in the values of ω that minimize the denominator. Now we are going to analyze the general expression for the denominator of the power spectrum, and show explicit expressions for the basal response.

The matrix A in the expression of $\Phi$ is the jacobian matrix which can be written in terms of the eingenvalues, then we can write:

$$\begin{eqnarray}\begin{array}{*{20}{l}}\text{det} \Phi(\omega)=\text{det}\left(\begin{array}{c} \text{i}\omega -{{\lambda}_{1}} & 0\\ 0 & \text{i}\omega -{{\lambda}_{2}} \end{array}\right)\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,=-{{\omega}^{2}}-\text{i}\omega \left({{\lambda}_{1}}+{{\lambda}_{2}}\right)+{{\lambda}_{1}}{{\lambda}_{2}}\end{array}\end{eqnarray} \tag{ A.18 }$$

$$\begin{eqnarray}&&\begin{array}{*{20}{r}}{{P}_{i}}(\omega)\propto \frac{1}{\left(-{{\omega}^{2}}-\text{i}\omega \left({{\lambda}_{1}}+{{\lambda}_{2}}\right)+{{\lambda}_{1}}{{\lambda}_{2}}\right)}\\ =\frac{1}{\left(-{{\omega}^{2}}+\text{i}\omega \left({{\lambda}_{1}}+{{\lambda}_{2}}\right)+{{\lambda}_{1}}{{\lambda}_{2}}\right)}\end{array}\end{eqnarray} \tag{ A.19 }$$

$$\begin{eqnarray}&&\begin{array}{*{20}{r}}{{P}_{i}}(\omega)\propto \frac{1}{{{\left({{\omega}^{2}}-{{\lambda}_{1}}{{\lambda}_{2}}\right)}^{2}}+{{\omega}^{2}}{{\left({{\lambda}_{1}}+{{\lambda}_{2}}\right)}^{2}}}\\ =\frac{1}{{{\left({{\omega}^{2}}-\text{det}(A)\right)}^{2}}+{{\omega}^{2}}{{\left(\text{Tr}(A)\right)}^{2}}}\end{array}\end{eqnarray} \tag{ A.20 }$$

This form for the power spectrum shows clearly the existence of a resonance: for a specific value of ${{\omega}^{2}}$ the denominator becomes small, and the power spectrum has a large peak centered on this frequency. The condition $\text{d}{{P}_{i}}(\omega)/\text{d}\omega =0$ gives:

$$\begin{eqnarray}&&2\left({{\omega}^{2}}-\text{det}(A)\right)2\omega +2\omega {{\left(\text{Tr}(A)\right)}^{2}}=0\end{eqnarray} \tag{ A.21 }$$

Then the condition (A.21) simply becomes: ${{\omega}^{2}}=\text{det}(A)-\frac{{{\left(\text{Tr}(A)\right)}^{2}}}{2}$, it implies that: $2\text{det}(A)>{{\left(\text{Tr}(A)\right)}^{2}}$, this condition in terms of the stability matrix implies that the eigenvalues of A are complex. It means that if the power spectrum has an extreme the stability matrix has complex eigenvalues.

Replacing in the expression for the power spectrum the eigenvalues for the model of the system in basal conditions: ${{\lambda}_{1}}=\frac{-{{k}_{2}}}{2}+\text{i}\sqrt{{{k}_{2}}{{k}_{4}}-\frac{k_{4}^{2}}{4}}$ and ${{\lambda}_{2}}=\frac{-{{k}_{2}}}{2}-\text{i}\sqrt{{{k}_{4}}{{k}_{2}}-\frac{k_{4}^{2}}{4}}$

$$\begin{eqnarray}&&{{P}_{i}}(\omega)\propto \frac{1}{{{\omega}^{4}}+{{\omega}^{2}}\left(k_{4}^{2}-2{{k}_{4}}{{k}_{2}}\right)+k_{4}^{2}k_{2}^{2}}\end{eqnarray} \tag{ A.22 }$$

$$\begin{eqnarray}&&P_{i}^{\prime}(\omega)\propto \frac{4{{\omega}^{3}}+2\omega \left(k_{4}^{2}-2{{k}_{4}}{{k}_{2}}\right)}{{{\left({{\omega}^{4}}+{{\omega}^{2}}\left(k_{4}^{2}-2{{k}_{4}}{{k}_{2}}\right)+k_{4}^{2}k_{2}^{2}\right)}^{2}}}=0\end{eqnarray} \tag{ A.23 }$$

${{P}_{i}}(\omega)$ has an extreme in $\omega =\sqrt{{{k}_{4}}{{k}_{2}}-\frac{k_{4}^{2}}{2}}$

As we can see, the imaginary part of the eigenvalues $\left(\text{Im}\left[{{\lambda}_{1,2}}\right]\right)$ and the value of frequency where ${{P}_{i}}(\omega)$ has an extreme are similar. Now we are going to write it in a more clearly way:

$$\begin{eqnarray}&&\begin{array}{*{20}{r}}\text{Im}\left[{{\lambda}_{1,2}}\right]=\sqrt{{{k}_{2}}{{k}_{4}}-\frac{k_{4}^{2}}{4}}=\sqrt{{{k}_{2}}{{k}_{4}}}\sqrt{1-\frac{k_{4}^{2}}{4{{k}_{2}}{{k}_{4}}}}\\ \approx \sqrt{{{k}_{2}}{{k}_{4}}}\left(1-\frac{k_{4}^{2}}{8{{k}_{2}}{{k}_{4}}}\right)\end{array}\end{eqnarray} \tag{ A.24 }$$

The above Taylor expansion may be a good approximation if we take into account that for the corresponding values of the parameters $\frac{k_{4}^{2}}{4{{k}_{2}}{{k}_{4}}}<\frac{k_{4}^{2}}{2{{k}_{2}}{{k}_{4}}}<0.5$.

In the same way:

$$\begin{eqnarray}&&\begin{array}{*{20}{r}}\omega =\sqrt{{{k}_{4}}{{k}_{2}}-\frac{k_{4}^{2}}{2}}=\sqrt{{{k}_{2}}{{k}_{4}}}\sqrt{1-\frac{k_{4}^{2}}{2{{k}_{2}}{{k}_{4}}}}\\ \approx \sqrt{{{k}_{2}}{{k}_{4}}}\left(1-\frac{k_{4}^{2}}{4{{k}_{2}}{{k}_{4}}}\right)\end{array}\end{eqnarray} \tag{ A.25 }$$

Summarizing, the power spectrum of the fluctuations has a maximum at a frequency close to the imaginary part of the eigenvalue that characterizes the frequency of the damped oscillations in the deterministic regime.