Noise Amplification in Human Tumor Suppression following Gamma Irradiation

The influence of noise on oscillatory motion is a subject of permanent interest, both for fundamental and practical reasons. Cells respond properly to external stimuli by using noisy systems. We have clarified the effect of intrinsic noise on the dynamics in the human cancer cells following gamma irradiation. It is shown that the large amplification and increasing mutual information with delay are due to coherence resonance. Furthermore, frequency domain analysis is used to study the mechanisms.

Mathematical models have achieved oscillatory dynamics by introducing ad hoc time delays to reproduce those that a system incurs when the various molecular components are manufactured [21,22,24,25,28,30,[34][35][36][37][38][39]. Related works have been performed on many fields of research, where delays were found to play a central role. For example, the importance of delay has also recently been recognized in neuronal dynamics [40][41][42][43]. From the mathematical point of view, the difference between single-cell experiments and cell population experiments of simple regulatory networks arises from stochastic events in individual cells that are averaged out in cell population. As the noise intensity of the regulating species increases, the noise intensity of the regulated one also appears to increase. Noise can induce many phenomena in nonlinear dynamical systems, including stochastic resonance, coherence resonance, pattern formation and so on. Lots of original research [44][45][46][47][48][49] and review [50][51][52] articles have been devoted to the stochastic resonance phenomenon. Noise-induced patterns in semiconductor nanostructures have been recently investigated by means of theoretical models [53], where random fluctuations play an essential role. Our presented results are crucially relying on coherence resonance, which has been recently studied for temporal systems [54][55][56][57] and spatially extended systems [58][59][60][61][62][63]. Specifically the relevance of intrinsic noise was elaborated on periodic calcium waves in coupled cells [64] and spatial coherence resonance in excitable biochemical media [65] induced by internal noise. A recent comprehensive review [66] has been done on the stochastic coherence. The large amplification results from the existence of coherence resonance with delay and noise.
In this article, by exploiting a microscopical signal-response model which was proposed in our previous articles [37,38] for studying the dynamical mechanism of the oscillatory behaviors for the activities of p53 and Mdm2 proteins in individual cells, we will explore the mechanism of noise amplification by considering the stochastic events in the cells.

Noise amplification
We introduce the probability Pr(n P ,n M ,t) for the p53 and Mdm2 populations P(t),M(t) ð Þ n P ,n M ð Þ. Then the master equation for Pr(n P ,n M ,t) is given by dPr(n P ,n M ,t) dt~P (n P ,n M )Pr(n P ,n M ,t) where t is added to account for the time delay between the activation of p53 and the induction of Mdm2.
Pr(n P ,n M ,t; m P ,m M ,t{t) is the joint probability distribution of having n P p53 molecules, n M Mdm2 molecules at time t and m P p53 molecules, m M Mdm2 molecules at time t{t. E P and E M are the unitary shift operators, E P Pr(n P ,n M ,t)~Pr(n P z1,n M ,t), E M Pr(n P ,n M ,t)~Pr(n P ,n M z1,t), and S P , a P , c P , m P , S M , a M , m M , K, N and S(t) are the parameters denoting various mechanisms as represented in our previous papers [37,38]. Assume that the time delay t compared with other characteristic times of the system is large, so the processes at time t and t{t are weakly correlated as Pr(n P ,n M ,t; m P ,t{t)~Pr(n P ,n M ,t) Pr(m P ,t{t). Adopting this approximation, we get dPr(n P ,n M ,t) dt~P (n P ,n M )Pr(n P ,n M ,t) The generating function G(s 1 ,s 2 ,t) is defined as G(s 1 ,s 2 ,t)~X We convert the infinite set of ordinary differential equations (3) to a single partial differential equation for G(s 1 ,s 2 ,t), LG The moments of the probability distribution can be found by expanding the generating function near s 1 ,s 2 ð Þ~1,1 ð Þ, LG LG n P ,n M~0 n P n M Pr(n P ,n M ,t)~SP(t)M(t)T, ð8Þ Substituting the expansion into Eq. (5) we obtain where the functions a 1 (t), a 2 (t) and b 12 (t) are Eqs. (6), (7) and (8), respectively. Above is the presentation of the derivation by help of generating functions. In fact, it delivers the same moment equations as the derivation by averaging the master equation. Both approaches run finally into equivalent approximations and problems if decoupling the moments. By the comparison between Eqs. (12) and the corresponding deterministic equations described in our previous papers [37,38], we find that due to the limit cycle of the deterministic description [37,38] changes to a decaying scheme as shown in Fig. 1.
From our numerical results, the reason for the decaying can be considered as dephasing that is mainly caused by differences in the Hill function P N (t{t)= K N zP N (t{t) ð Þbetween the cells. The reason that Hill functions are different is the different states of the different cells at time t{t, i.e., some dephasing happened at time t{t for it to have this impact. The delay further amplifies the differences between cells, causing further dephasing. but if we take two cells with identical state space paths, their Hill functions will also be the same.
This initial difference between the particle numbers of chemical species in different cells, which causes the difference in Hill function at later time, is entirely caused by the intrinsic noise. In fact, any oscillating chemical system, with or without delayed dynamics, will demonstrate dephasing between different realizations, and it isn't an artifact of the delayed dynamics themselves, although this will undoubtedly cause further decorrelation of different realizations at later time, which causes the damped behavior at the population level (which can be thought of as simply taking a large number of realizations of the same stochastic system). Essentially, the value which the cell population converges to is simply approximately the mean of the invariant distribution of the chemical species for one cell, multiplied by the number of cells in the population of interest. This can be shown more rigorously for large populations using the ergodic property of the system. Fig. 2 shows the average power spectrum S P (v) for P(t) time series as a function of frequency v. We also plot the spectrum of the corresponding deterministic model, with delay (e.g., time delay t d~1 00 min in Fig. 2) but without noise, to compare its spectrum with stochastic ones. It can be clearly seen that S P (v) without noise is much smaller than those with noise. Significantly, for the cases with large t (especially, t is larger than the Hopf bifurcation point t c ), there are obvious peaks appearing in S P (v) for P(t) at v=0. This tells us that there is a very large amplification of intrinsic noise due to the resonant effects. This characteristic phenomenon may be termed as coherence resonance with delay and noise, for distinguishing from the '' stochastic resonance'' in common sense.
The peak frequency corresponds to the characteristic frequency of the solution of Eqs. (12), which represents the mean frequency of Fourier transform F P(t) ½ . It is very intriguing that the width of S P (v) represents the dephasing effects, which gives the damping strength on the amplitude of vP(t)w. In order to analyze this resonant oscillation more transparently, we phenomenologically fit S P (v) for the cases with large t (twt c ) shown in Fig. 2

by a formula
where the parameters a, b, V and C are t-dependent. Note that Eq. (14) can be analytically derived with the chemical Langevin equations corresponding to Eqs. (1) under the linearization approximation.
The resultant V and C are shown in the inset picture of Fig. 2. It is obvious that the mean frequency V decreases against t, which is consistent with the conclusion described in our previous article [37]. This is particularly important in biology because in general the low frequency is much more significant than higher frequency in biological systems. C also decreases as t increasing, which means that the oscillation may dominate the evolution of vP(t)w and lasts for rather longer time for very large t. This phenomenon is very intriguing from the biological point of view because it may tell us that the time delay induced by the underlying multistage reactions may weaken the effects of stochasticity and strengthen the oscillation of the relevant molecules.
Mutual information (MI) is meaningful to discuss resonant phenomena [67], so we give the mutual information between the two components p53 and Mdm2 in the nonlinear delayedfeedback network motif. MI is a measure of the amount of information that one random variable interacts with another. It is the reduction in the uncertainty of one random variable due to the  MI is zero if and only if the two random variables are strictly independent [69]. Numerically calculating the mutual information between trajectories is in general a formidable task [70], since the joint distribution of continuous variable is smoothly obtained only for large scale stochastic simulation. Intensive work has been done on estimating the mutual information. Khan et al. [71] reviewed three MI estimators: Kernel density estimators, k-nearest neighbor method and Edgeworth expansion. Recently, Suzuki et al. [72] proposed a novel MI estimator called Least-Squares Mutual Information, and discussed the characteristics of the three existing approaches. However, it is accessible here due to the discreteness of the system with the exact delay stochastic simulation algorithm (DSSA) [73]. Information theory [74] provides a natural framework for many problems in biological information processing. The Shannon mutual information has been applied to study the stochastic resonance (SR) [67,75,76], instead of the signal-tonoise ratio (SNR). It can be seen from Fig. 3 that when the DNA is damaged, the phosphorylation of p53 modifies its binding properties to Mdm2, so MI is small; But when the signal is completely resolved, e.g., after t th~1 750 min, MI is large because the amount of p53 is kept low and tightly regulated by the genetic network of Mdm2 and p53 itself. Fig. 4 shows that MI in steady state increases with the increase of time delay due to the coherence resonance.

Fourier analysis
To describe the nonlinear dynamics more clearly, we use frequency domain analysis method to study the mechanisms of the p53 network motif. Our model can be described by a set of chemical Langevin equations corresponding to Eqs. (1), dP(t) dt~S P {a P M(t)P(t)(1{c P S(t)) where g 1 (t) and g 2 (t) are Gaussian white noise, Sg i (t)T~0, Sg i (t)g j (t 0 )T~Sg i g j Td(t{t 0 ), i,j~1,2 f g [77]. In order to analyze our model in the frequency domain, we first replace P(t) and M(t) in Eqs. (18) by where P Ã and M Ã represent the stationary solutions of the deterministic equations of Eqs. (18) with t~0, which satisfy the equations Since we are discussing the solution in the oscillatory scheme, here the signal S(t) is set to be 1. If one hopes to discuss the case of the  stationary solution in t??, he can simply set the parameter c P ?0 mathematically, because at t??, the damage can be supposed to be completely resolved as S(t??)~0, i.e., the signal S(t) is first set to 1, later c P is removed because S(t) is becoming 0 if time tends to infinity. Then Eqs. (18) can be rewritten as dp dt~A where and the nonlinear term is kept up to the second order in p(t{t). The Fourier transformations of Eqs. (18) take the form where Since Eqs. (23) are integral equations, they can be solved by interpolation method and truncated at a specific order, the following calculation includes convolutions in the spectral presentation replacing the nonlinear items in the temporal one and truncating them, e.g., we can first solve the linear equation substitute the solutions p(v) and m(v) of Eqs. (25) into Eqs. (24), and then F (v) and G(v) are functions of v. Under the approximations of weak noise and weak negative feedback mechanism, in this paper, the solutions of both p(v) and m(v) are retained up to the second order of g 1 (v) and g 2 (v), because for Gaussian noise, the terms of higher order can be omitted in Ito-Wiener approximation. The validation of such approximations will be discussed with our numerical simulation later. We define and then it can be derived from Eqs. (23) that By defining the intermediate variables, I 4 (v,v 0 )~g 2 (v 0 )g 4 (v 0 )zg 5 (v 0 )zg 2 (v 0 )g 6 (v) ð Þ g 2 (v{v 0 ),ð29dÞ Eqs. (28) can be written as {?
where k,m,n~1,2,3,4 f g . Then the correlation functions of p(v) and m(v) can be expressed as S m (v)~Sm(v)m Ã (v')T b 1 Sg 2 1 Tzb 2 Sg 1 g 2 Tzb 3 Sg 2 2 T zb 4 Sg 2 1 T 2 zb 5 Sg 2 2 T 2 zb 6 Sg 1 g 2 T 2 zb 7 Sg 2 1 TSg 1 g 2 Tzb 8 Sg 1 g 2 TSg 2 2 T: The parameters a 1 , b 1 , a 2 , b 2 , Á Á Á, a 8 , b 8 represent the contributions of Sg 2 1 T, Sg 1 g 2 T, Á Á Á, Sg 1 g 2 TSg 2 2 T to the correlation functions S p (v) and S m (v), respectively. With the aid of the intermediate variables, those parameters can be expressed as 3,4 (v)zJ 4,3 (v)zL 2,4 (v)zL 4,2 (v)zL 3,4 (v)zL 4,3 (v)Þ, For S p (v), For S m (v), It is worthwhile to mention that a module, which consists of two components, has been discussed recently [78]. They studied a set of coupled Langevin equations for the interacting species. It is very interesting that in the absence of delay and nonlinearity, i.e., a special case of the spectrum as t~0 in Eqs. (55), Eqs. (34) can be reduced as which are consistent with the results presented in the previous paper [78]. Another characteristic feature of Eqs. (34) is that when Sg 1 g 2 T is assumed to be zero, which means that g 1 and g 2 are uncorrelated, both S p (v) and S m (v) can be written as a sum of two contributions which is the so-called spectral addition rule as derived in the previous paper [78]. Even in this case, the coefficients in our results still include the effects coming from the time delay and negative feedback mechanism.
In our numerical calculation, we use the fourth-order stochastic Runge-Kutta method for integrating the chemical Langevin equations (18), and Gaussian integration method to calculate the integrations in Eqs. (34). The numerical results have shown that the correlation functions S p (v) and S m (v) for p(t) and m(t) are precisely consistent between the ones with chemical Langevin equations (18) and the ones with Eqs. (34), which verifies our truncation method in Eqs. (23). The Fourier transforms of p53 and Mdm2 dynamics show that the number of the resonant peaks would increase as time delay increases, which is consistent with the experimental results [12]. The general finding of our analysis is that an increase of delay between activation and induction induces an oscillatory behavior with frequency which corresponds nearly to the delay time. The spectral analysis as well as the mutual information supports this finding. The general finding is in good agreement with our previous work [37].
Bioscience and nanoscience provide pretty examples of nonequilibrium and nonlinear dynamics in which noise can be expected to have unavoidable effects. The methods developed over years to deal with the effects in physical systems will help us to further our understanding of the mechanisms ascribed to nonlinearity and noise.

Methods
The stochastic p53 circuit was characterized by a Monte Carlo method called the exact DSSA. Numerical integration of the equations was carried out using Matlab software.