Period doubling induced by thermal noise amplification in genetic circuits

Rhythms of life are dictated by oscillations, which take place in a wide rage of biological scales. In bacteria, for example, oscillations have been proven to control many fundamental processes, ranging from gene expression to cell divisions. In genetic circuits, oscillations originate from elemental block such as autorepressors and toggle switches, which produce robust and noise-free cycles with well defined frequency. In some circumstances, the oscillation period of biological functions may double, thus generating bistable behaviors whose ultimate origin is at the basis of intense investigations. Motivated by brain studies, we here study an “elemental” genetic circuit, where a simple nonlinear process interacts with a noisy environment. In the proposed system, nonlinearity naturally arises from the mechanism of cooperative stability, which regulates the concentration of a protein produced during a transcription process. In this elemental model, bistability results from the coherent amplification of environmental fluctuations due to a stochastic resonance of nonlinear origin. This suggests that the period doubling observed in many biological functions might result from the intrinsic interplay between nonlinearity and thermal noise.

T he diversity of multicellular organisms is regulated by a system of rhythms that affect all forms of their functionalities [1][2][3] . A feature shared by many organic systems, however, is the period doubling of their internal cycles. Appearing spontaneously or in response to external disturbances, the cycle of cellular organisms is observed to become more complex, enriching its activity and leading the system to complete its functionality with a period that is two times larger than initial one [4][5][6][7][8][9][10][11] . This phenomenon is observed in a wide range of scales and can be found, for example, in the arrhythmia of cardiac functions 4 , or in the firing of neurons 5,6 , where it has been tentatively explained on the basis of the well know bifurcation to chaos paradigm 9,10 , or in the distribution of population growth in yeast 12 , where it was ascribed to genetic regulatory circuits. Period doubling has also been observed during cell divisions and circadian cycles, where it can manifest spontaneously 11 or in response to chemical perturbations 13 . Similarly, it has been shown that the genome-wide transcriptional oscillation in yeast can experience period doubling in reaction to drugs 14 . The ubiquity of this phenomenon challenges the stability of the ''biological oscillators'' 15,16 , and still deserves a comprehensive and clear interpretation.
In this paper, we investigate the doubling process starting from an elementary genetic circuit that combines two ubiquitous aspects of biosystems: inhibitory mechanisms and noise. The active role of noise in genetic and biochemical circuits has already been emphasized [17][18][19][20][21][22] and its amplification has been recently found at the basis of novel and important physical phenomena 23 . On the other hand, the importance of inhibitory functions has already been recognized in neurobiology, where they originate nonlinear networks that are at the basis of a diversity of brain functionalities [24][25][26] . We demonstrate that our elemental genetic circuit can be mapped on an electric RLCD series, made by a resistor, an inductance, a capacitor and a diode, which represent a general model of resonant wave interactions in the presence of a strongly nonlinear response. In this elementary model, nonlinearity arises from the inhibitory function of the diode, which regulates its voltage (i.e., the protein concentration in the genetic counterpart). We show that the two aspects of the genetic dynamics, namely nonlinearity and noise, do not produce any interesting evolution if considered alone. On the contrary, when noise and nonlinearity are combined together, period doubling appears. Analytic theory demonstrates that such phenomenon is sustained from nonlinear parametric resonance, which strongly amplifies noise in a narrow band that resonates with the coherent part of environmental fluctuations. This suggest that period doubling commonly observed in biological activity can be ascribed to the interplay between nonlinearity and noise.

Results
From genetic models to nonlinear electric circuits. We consider a simple protein transcription process described by the functional autorepressor 27 sketched in Fig. 1a, with its associated dynamics for the protein p(t) and mRNA m(t) concentrations, respectively 28,29 : In Eqs. (1), a i and b i stand for the reaction and degradation rates, respectively, while g(p) is the promoter-inhibitor activity function. The latter models the mechanism where the gene product represses its own transcription. The function g is an upper bounded decreasing function of p 30 , which can be expanded in series: with p 0 a constant parameter. We begin our discussion by considering the simplest case given by a first order function g(p)5(p 0 2 p), and then discuss the effects of high order terms (see the Discussion Section). Reaction/degration rates are complemented by the presence of an environment at nonzero temperature, modeled by the time dependent terms S p [p(t)] and S m (t). The function S p , in particular, accounts for an additional protein concentration depletion due to its consumption in the different cell compartments, while S m (t) models the action of the environment on the mRNA transcription rate 2 . In order to pursue a general theory, we decompose S m into two main contributions S m (t) 5 s(t) 1 g(t). The term s(t) represents the coherent (periodic) contribution originating from cyclic mechanisms (such as e.g. circadian rhythms and/or cell divisions), while g(t) models the incoherent (noise) part arising from random fluctuations. The fluctuations arising from S m , due to the coupling between m and p, provide a statistical noise source also for the protein concentration p. We finally model the term S p by a nonlinear p dependent degradation function, which we express in the simple form: being x p the nonlinear degradation rate of the protein concentration. The nonlinear response described by Eq. (3) is illustrated in Fig. 1c (solid line) and models in a simple fashion the mechanism of cooperative stability, which is widely observed in many experiments 31,32 and leads to an enhanced degradation rate at low protein concentration and a lower depletion at higher concentration 33 . Equation (3) has the dynamics of a ''biological'' diode, which leads to an interesting equivalent electric circuit representation. In order to illustrate the circuit analogy within the simplest theoretical framework, we model the exponential response in Eq. (3) by a piecewise linear model ( Fig. 1c dashed line). This type of approximation is largely employed in electronic circuits to provide a simple yet accurate representation of nonlinear components such as diodes and transistors (see e.g., chapter 3 of 34 ). The piecewise linear model is characterized as follows: with a constant d p estimated from (3): d p~xp logpz1 ð Þ, beingp the characteristic scale for the protein concentration in the autoregressor dynamics. For the situation illustrated in Fig. 1c with x p 5 0.2, if we assume a characteristic concentrationp<100 we obtain d p 5 0.9. Equations (4) generate the following equations of motion when the protein concentration reaches p 5 0 (OFF state): where only the mRNA m(t) is allowed to evolve while the protein p concentration stays in its minimum (zero) value. The duration of this state depends on the value of the time derivative _ p in Eq. (1). When the latter is negative, and therefore showing a tendency of the system to decrease the protein concentration below the minimum, the dynamics is maintained in the OFF condition. When _ p becomes positive, conversely, the system gets back to the ON state and the concentration p increases again. From the second of Eqs. (1), we immediately observe that the threshold for a positive protein production, when p50, is represented by a p m(t) 5 S p (i.e., when a p m(t) $ S p we have _ p t ð Þ §0). It is worthwhile emphasizing that such ON-OFF states provide a piecewise linear representation for Eq.
(3) and do not generate any discontinuity in the bio-physical system, which follows the continuous dynamics described by Eqs.
with V f and C the forward bias and the junction capacity of the diode, respectively. The circuit constants (R, L, C, V f , V(t)) and the  dimensionality scaling constant c 0 are related to the autorepressor parameters by: Under the coordinates change (6) -(7), the functional autorepressor (1) -(5) is mapped into the equivalent electric dynamics: when the diode is conducting (ON state). Switching conditions are promptly interpreted from (6) as the conduction/not conduction states of the diode, and read OFFRON when y . 0 and ONROFF when x5CV f . These are equivalent to the conditions a p m . d p and p50 of the genetic circuit (1) - (5). To account for the full nonlinear form of the protein depletion function in Eq. (3), we need to consider the exponential nonlinear response of the diode in the circuit of Fig. 1c, thus obtaining: : being y s the current flowing into the diode and V T its equivalent thermal voltage. The latter can be estimated by inverting Eq. (11), obtaining V T~Vf logỹ s z1 À Á , withỹ s the characteristic current flowing into the diode. For a typical current y s < 1mA and a potential V f 5 jxj/C < 0.6 V, we have V T < 600 V. The protein nonlinear depletion rate x p characterizes the equivalent thermal voltage V T of the diode through Eqs. (7) that, in turn, define the diode potential V f . Equations (10) -(11) constitute the equivalent electric circuit of Eqs. (1) and (3), with Eq. (11) being the electric counterpart of (3). Equations (5) and (9), conversely, represent the corresponding piecewise linear models. Given a set of ''genetic'' parameters a i , b i , S i , the equivalent electric system possesses two free constants: L and V f . While the diode bias V f defines the current amplitude scale, the inductance L provides an arbitrary scale for the resonant frequency v 0~1 . ffiffiffiffiffiffi LC p and quality factor Q 5 v 0 L/R of the circuit when the diode is conducting. Besides the present context where the RLCD circuit has been derived, we observe that Eqs. (10) -(11) yield a general model of resonant wave interactions in the presence of a strongly nonlinear response [i.e., Eq. (10)] and a generic excitation V(t) constituted by both periodic and random terms. This model can be interesting in many different contexts and in particular in nonlinear optics, where oscillators similar to (10) - (11) are at the base of the study of various types of light-matter interactions 35,36 . This aspect will be systematically investigated in future work. In order to highlight the simple nature of the coordinates transform (6) - (7), we cast Eqs. (1) and (9) in the form _ x~J : x~S, with x 5 (m, p) and: for the biological case, and x 5 (x, y) with: for the electric case. In absence of external sources [i.e., for S m (t)50 or V(t)50], Eqs. (12) -(13) possess an equilibrium point, which we indicate by ( m, p) or ( x, y). Such a fixed point is promptly calculated from _ x~0 and reads The simplicity of the coordinate transformation (6) - (7) is now clear by comparing the two expressions in Eqs. (14), which show that the RLCD system offers a simple representation where the fixed point lies in the origin of the coordinates. The scaling constant c 0 can be recast in the form c 0~C V f p, thus indicating that the divergence of c 0 for a vanishing denominator is approached when p?0 from above, i.e. when the equilibrium point is at the limit of the field of existence of p.
Numerical results with sinusoidal inputs. We studied the RLCD system via numerical integration of their respective equations of motion. Figure 2 summarizes our results in the case of a fluctuating sinusoidal periodic source V(t)5V 0 sin(v p t) 1 g(t), being g a gaussian noise term with zero average AEgae50 and time uncorrelated (white) variance AEg(t)g(0)ae52 T L d(t). The noise variance is measured by the ''temperature'' T, which represents the total spectral power of the fluctuations. In our simulations, we investigated the effect of incoherent (Fig. 2a-c) and coherent ( Fig. 2d-f) contributions separately, as well as their combination (Fig. 2g-i). In order to simulate a realistic set of electric parameters, we considered L 5 0.01 H, C 5 310 213 F, R52 V, V f 5 20.6 V and jV T j 5 600 V. This choice gives a resonant frequency v 0 5 0.58 MHz and a quality factor Q 5 2887. In the example reported here, the source is characterized by V 0 5 1.6 V, v p 5 1.4v 0 , and T 5 10 27 W. Figure 2a-c displays the system dynamics when only incoherent noise g(t) is present, with temperature T 5 10 27 W. At this temperature, random fluctuations have small amplitudes, and are not strong enough to generate any current in the diode. The circuit therefore behaves as a simple RLC series, with bare frequency v 0 and damping v 0 /Q, and its main effect is to produces a ''colored'' noise ( Fig. 2b-c). In the genetic counterpart, this implies that the system maintain a simple oscillating protein concentration p(t). When only a coherent periodic signal V o sin(vt) is considered (Fig. 2d-f), with the present parameter choice the source intensity is sufficiently large to induce current in the diode. In this regime, circuits similar to that depicted in Fig. (5) have been investigated in the past, both experimentally [37][38][39][40] and numerically 41,42 . It has been found that when the diode is characterized by a fixed junction capacity C and when there is no delay in the diode switching -as in the case studied here-the system displays a simple resonant behavior, without bistability, bifurcation or chaos. This behavior is confirmed in the present study: as shown in Fig. 2f, the power density spectrum P v ð Þ~ỹ v ð Þ j j 2 (hereỹ v ð Þ is the Fourier transform of y(t)) is particularly simple, it only shows a strong peak at the input frequency v 5 v p and a rectified component at v 5 0, which originates from the rectifying action of the diode. The very interesting dynamics is observed in Fig. 2g-i, when both coherent and incoherent sources are simultaneously present. In this www.nature.com/scientificreports case, a strong harmonic component at a new frequency, v 5 v p /2, appears, which represents a subharmonic of the input frequency v p , and period doubling is observed in the dynamics (Fig. 2h-i). Quite remarkably, the strong noise resonance at v 5 v 0 (Fig. 2c) disappeared from the dynamics. By comparing the spectra in Fig. 2i and Fig. 2f, we observe that the current power densities P(v) at v 5 v p (i.e., the input pump frequency) and v 5 0 (i.e., the rectified component) are approximatively the same. The strong subharmonic peak at v p /2 is therefore the result of the amplification of the noise fluctuating in the background. Stochastic contributions, in particular, are coherently amplified and constructively interact with the input signal thereby sustaining a period doubling (Fig. 2h). To further investigate this process, we calculated the intensity of the peak at v 5 v p /2 for different input frequencies v p and input amplitudes V 0 , as reported in Fig. 3. The subharmonic peak at v p /2 has a resonant-like intensity as a function of the pump frequency, having its maximum when v p < 3v 0 /2 for low value of V 0 and red-shifting at higher pump intensity. On increasing V 0 , the subharmonic peak maximum increases, disappearing for V 0 below a v p -dependent threshold. From Fig. 3, we observe that the efficiency of the subharmonic generation process strongly increases in the region where v p < v 0 .
Period doubling in the presence of short pulses. Periodic rhythms of biosystems are often observed as sequences of spikes, with each spike characterized by a short-living pulse. Is therefore important to investigate the occurrence of period doubling in this specific case. Figure 4a-c summarizes our results by considering square pulses of time duration dt/T 5 1/6, with frequency v p 5 2p/T and amplitude V 0 5 2.2 V. The noise temperature is set to T 5 5 ? 10 27 W. Despite the different time evolution of the source, period doubling is qualitatively identical to the one observed for sinusoidal inputs, thus witnessing the robustness and ubiquity of the phenomenon. From a quantitative perspective, the only appreciable difference is observed in the spectrum (Fig. 4c), where we report a larger bandwidth of the peak at v 5 v p /2 with respect to the sinusoidal case (Fig. 2i). We conclude our numerical campaign by investigating the period doubling dynamics versus the noise temperature. To this extent, we calculated the current y for different temperatures T, and we extracted the current minima y n 5 y(np/v p ), which are obtained by sampling the current at half of the input period p/v p . For each temperature T, we averaged over 40 realization of noise and collected many sequences of y n , in order to obtain a significantly large statistic, which we reported in Fig. 4d. The figure shows that the ''stochastic'' period doubling observed in Fig. 2-4 does not belong to the classical bifurcation scenario of chaotic systems: no abrupt change in the period takes place, and the development of the doubling is a gentle process without any clear threshold in the control parameter (the temperature T). As seen from Fig. 4d, the distance between the current minima increases nonlinearly with the temperature, witnessing the strongly nonlinear nature of the process.  Discussion In order to understand the role of the non linearity and, in particular, the physics of the period doubling observed in Fig. 2-4, we consider the piecewise linear model (9) and refer to Fig. 5, which displays the orbit of the system in the (x, y) plane for sinusoidal sources (blue line) and short pulses (green line) with and without noise. In absence of random fluctuations (Fig. 5a), the nonlinearity of the dynamical evolution is clearly seen as a discontinuity in the orbit tangent, taking place at points A, E (the point of the ONROFF shift) and B (OFFRON shift). When the orbit reaches points A or E the charge x is kept equal to CV f until the current y become newly positive in B.
Following the diode in the ON state, for a sinusoidal source the orbit moves first along the semicircle (B-D), and then along the quarter of circle (C-A), which completes the cycle. In the case of square pulses, the dynamics is the same with the only difference of a smaller cycle B-D-E-B. The presence of noise does not change qualitatively the picture, but spreads out the orbits in a random fashion in both cases (Fig. 5b). Is worthwhile observing that part of the plane x , CV f , even for y . 0, is never visited (at point A or E, in fact, the charge variation _ x is always positive). This observation allow us to write Eqs. (9) in a more convenient form, where the switching ON-OFF conditions involve only the variable x: We can now study Eqs. (15) to investigate the origin of period doubling. To this aim, we employed a perturbation analysis from the solution of the nonlinear system in the absence of noise, i.e., for T 5 g 5 0. Due to the smallness of the subharmonic peak with respect to the input frequency at v p (Fig. 2i, 4c), we apply first order perturbation theory. In particular, we indicate with [j(t), f(t)] the nonlinear, noise-free solution of Eqs. (15) for g(t) 5 0. Although we cannot write this solution in closed form, we know that this is strongly oscillating at v 5 v p in time, and lies in the curve A-B-C-A (or B-D-E-B) of Fig. 5a. In the presence of nonzero noise, we set x(t)5j(t) 1 D x (t) and y(t)5f(t) 1 D y (t). By substituting the latter expressions into Eq. (15), and retaining only the first order terms in D x (t) and D y (t), we obtain the following dynamics: where: being d(x) the Dirac-d function.
During an orbit in the (x, y) plane, during the ON state we can distinguish between two different regions for the noise-free solution j(t). The first region is characterized by j 1 D x . x o , which is represented by the segment B-C-A (or B-D-E) in Fig. 5a, with points  B, A (or B, E) and a small interval around them excluded. In this case h 5 1 and f 5 x, which yields h9 5 0 and f 9 5 1. Equations (16) can be trivially solved in the frequency domain v and the solution for the 2D y e ivt zc:c: (c.c. stands for complex conjugate) is a damped harmonic oscillator, with maximum amplitude at v5v 0 beingg the Fourier transform of g(t). The second region of the ON state, which is the most critical, is in the vicinity of point A and is represented by j{x o~, with =D x . In this situation, the term D x leads the dynamics to continuously oscillate between the A-B (or E-B) segment (ON state) in Fig. 5a and the rest of the cycle B-C-A (or B-D-E) (OFF state). When the system switches between these two segments, the charge D x experiences a discontinuous dynamics originated from the term h'ð Þ~2dð Þ, with rapidly changing from zero to a small but finite value due to the noise fluctuations. The current D y , conversely, evolves smoothly thanks to f 'ð Þ~1. In this condition, equations of motion (16) become: The term f(t) appearing in the RHS of the first equation oscillates with frequency v p , and this imposes a resonance condition with the oscillation frequency of D x (t). In the Fourier domain, in fact, the only frequency admissible is the one where the term f t ð ÞD x1 4D x e ivt zc:c: f e ivpt zc:c: oscillates at the same frequency with D x and D y . This condition imposes phase-matching v p 6 v 5 v, which is equivalent to v 5 v p /2. The only oscillation that can be observed in the evolution of D x , D y is therefore at v 5 v p /2, cause all the frequencies that do not satisfy phase-matching are not allowed in the second region and ruled out from the dynamics. Equations (19) allow to qualitatively assess the generality of our findings with respect to the particular shape of the promoter activity function g(p). In the most general situation, the function g(p) can be expanded in Taylor series following Eq. (2). Linear terms / p 0 and / p give rise, in the electric circuit representation, to linear dissipative electric components RLC. High order nonlinear terms / p n with n $ 2, conversely, originate nonlinear dissipative terms that do not qualitatively affect the existence of the parametric resonance at v p /2. The latter, in fact, originates from the nonlinear response of the diode that, as shown in Eqs. (6) - (7), is sustained by the mechanism of cooperative stability and the linear first order terms in the expansion of g(p). In order to solve for the currentD y , we represent d(x) 5 H ? rect(Hx), being rect(x) 5 h(x 1 1 2 2) 2 h(x 2 1/2). After taking the limit for H R ', we obtain the following solution: The spectral power J y~Dy v p 2 À Á 2 of the current at v 5 v p /2 can be then obtained by combining in time the two expressions forD y appearing in Eqs. (20) and (18) for v 5 v p /2 with weights 1 2 a and a respectively. The latter indicates the time spent in the different region of the cycle and critically depends on the dynamics of the noise free solution (j, f). When we combine in time different spectra with different time duration, the bandwidth and the Q-factor of the single spectra change as well, as we are convolving the spectrum arising from an infinite signal with the Fourier transform of a box function of a finite length. We can therefore define an effective Q-factorQ, which here must be consider as a fitting parameter as it depends on the time spent by j and g in the different regions of the cycle. After some straightforward algebra, we obtain the following expression for the current power density: being V 5 v p /2v 0 , jg r /Rj 2 the current noise spectral power at the peak frequency v p /2. Equation (21) has a nonlinear stochastic resonance at V < 0.5, which predicts a strong noise amplification of the subharmonic v p /2 when v p < v 0 . In order to verify our theory, we calculated the behavior of J y versus V by a series of simulations with a sinusoidal input with a varying frequency v p . For each v p , we averaged over 40 realization to calculate mean value and standard deviation of J y . We then compared numerical simulations versus Eq. (21), with parameters a 5 0.51 andQ~3:4 given by a nonlinear least square fit (Fig. 6a). We observed an excellent agreement between Eq. (21) and the results from numerical simulations, confirming the strongly resonant nature of the process. Equation (20) allows also to interpret the quantitative differences between the power density spectra in Fig. 2i and Fig. 4c. According to the phase matching condition imposed by Eq. (20), in fact, the bandwidth of the amplified noise is expect to match the bandwidth of coherent part the input source near v p , where each component resonate with its subharmonics and get amplified through parametric resonance. This process is highlighted in Fig.6b, where we superimpose different part of the current power density spectrum (Fig. 6b) obtained from a numerical simulation with V 0 5 1.6 V, T 5 10 26 W and v p /v 0 5 1.4. As seen in the figure, the power density of the amplified noise near v p /2 ( Fig. 6b solid green line) matches very well the spectrum of the input signal near v p (Fig. 6b solid red line). In the case of short pulses, parametric resonance is triggered by a larger spectrum near v p , due to the larger harmonic content of a short pulse with respect to a purely sinusoidal source, and therefore results into a larger amplified band near v p /2.

Beyond gene expression
The ''stochastic'' parametric amplification reported in this work provides a new panorama in nonlinear dynamics, whose importance goes also beyond the biological circuits introduced in this paper. A deep discussion on this subject goes clearly beyond the scope of this article, however a few interesting points can be highlighted. The nonlinear dynamics of a single non-chaotic oscillator possessing a resonant frequency at v 0 , is normally observed as the generation of different higher harmonics 2v 0 , 3v 0 , 4v 0 , … 43  generation of v 0 /2, or equivalently period doubling, is impossible as it does not follow from any integer combination of higher harmonics. Period doubling is typically observed, in fact, in chaotic resonators and is a well-known mechanism for routing the dynamics to chaos 44 . Our works open a new interesting scenario, where period doubling can be triggered by the noise interacting with a single non-chaotic oscillator. This opens the possibility to develop applications where noise acts as an active pathway for nonlinear dynamics, enabling functionalities that would otherwise be impossible to achieve for the system. Among the possible systems that might benefit from this process, we here mention the important case provided by midinfrared energy harvesting, where light-matter interactions are modeled by the same circuit of Fig. 1b 45 .
Conclusions. In conclusion, we have investigated the dynamics of a elemental genetic circuit characterized by the interplay between nonlinearity and noise arising from environmental fluctuations. We showed that the circuit can be mapped into an electric analogue represented by a RLC series with a diode D. In this elementary system, nonlinearity arises naturally from the diode, which acts as a inhibitory feedback that maintains the dynamics in equilibrium. When noise and nonlinearity are considered independently in the system, we do not observe any interesting evolution, while when they exist together we observe period doubling triggered by parametric amplification of thermal noise. Parametric resonance, in particular, amplifies noise in a specific band that resonates with the coherent part of the fluctuations of the environment. This result may explain the ubiquitous period doubling phenomenon reported in different cases in biological matter, and suggests that it may be traced back to the interaction between nonlinearity and environmental fluctuations. This works also highlights the active role of noise in inhibitory processes and suggests that its modeling is key to understand the complexity of biological functions, ranging from brain activity to gene expressions.