The diversity of multicellular organisms is regulated by a system of rhythms that affect all forms of their functionalities[@b1][@b2][@b3].  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[@b4][@b5][@b6][@b7][@b8][@b9][@b10][@b11].  This phenomenon is observed in a wide range of scales and can be found, for example, in the arrhythmia of cardiac functions[@b4], or in the firing of neurons[@b5][@b6], where it has been tentatively explained on the basis of the well know bifurcation to chaos paradigm[@b9][@b10], or in the distribution of population growth in yeast[@b12], 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[@b11] or in response to chemical perturbations[@b13].  Similarly, it has been shown that the genome-wide transcriptional oscillation in yeast can experience period doubling in reaction to drugs[@b14].  The ubiquity of this phenomenon challenges the stability of the "biological oscillators"[@b15][@b16], 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[@b17][@b18][@b19][@b20][@b21][@b22] and its amplification has been recently found at the basis of novel and important physical phenomena[@b23].  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[@b24][@b25][@b26].  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[@b27] sketched in [Fig.  1a](#f1){ref-type="fig"}, with its associated dynamics for the protein *p*(*t*) and mRNA *m*(*t*) concentrations, respectively[@b28][@b29]: In Eqs. ( 1), *α~i~* and *β~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*[@b30], 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*) = (*p*~0~ − *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[@b2].  In order to pursue a general theory, we decompose *S~m~* into two main contributions *S~m~*(*t*) = *σ*(*t*) + *η*(*t*).  The term *σ*(*t*) represents the coherent (periodic) contribution originating from cyclic mechanisms (such as e.g. circadian rhythms and/or cell divisions), while *η*(*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 *χ~p~* the nonlinear degradation rate of the protein concentration.  The nonlinear response described by Eq. ( 3) is illustrated in [Fig.  1c](#f1){ref-type="fig"} (solid line) and models in a simple fashion the mechanism of cooperative stability, which is widely observed in many experiments[@b31][@b32] and leads to an enhanced degradation rate at low protein concentration and a lower depletion at higher concentration[@b33]. [ Equation (3)](#m3){ref-type="disp-formula"} 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](#f1){ref-type="fig"} 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[@b34]).  The piecewise linear model is characterized as follows: with a constant *δ~p~* estimated from (3): , being the characteristic scale for the protein concentration in the autoregressor dynamics.  For the situation illustrated in [Fig.  1c](#f1){ref-type="fig"} with *χ~p~* = 0.2, if we assume a characteristic concentration we obtain *δ~p~* = 0.9. [ Equations (4)](#m4){ref-type="disp-formula"} generate the following equations of motion when the protein concentration reaches *p* = 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 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 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 *p* = 0, is represented by *α~p~m*(*t*) = *S~p~* (i.e., when *α~p~m*(*t*) ≥ *S~p~* we have .  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. ( 1) -- (3) at every instant.  To map Eqs. ( 1) -- (5) into a simple nonlinear electric circuit, we introduce the following linear change of variables (*p*, *m*) → (*x*, *y*), being (*x*, *y*) charge and current, respectively, of the RLCD circuit of [Fig.  1b](#f1){ref-type="fig"} with: 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 *γ*~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: with *V~d~*(*t*) = *V~f~* when the diode is in the OFF state and 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 OFF→ON when *y* \> 0 and ON→OFF when *x* = *CV~f~*.  These are equivalent to the conditions *α~p~m* \> *δ~p~* and *p* = 0 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](#f1){ref-type="fig"}, thus obtaining: with: 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 , with the characteristic current flowing into the diode.  For a typical current *y~s~* ≈ 1*mA* and a potential *V~f~* = \|*x*\|/*C* ≈ 0.6 *V*, we have *V~T~* ≈ 600 *V*.  The protein nonlinear depletion rate *χ~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)](#m10){ref-type="disp-formula"} -- [(11)](#m11){ref-type="disp-formula"} constitute the equivalent electric circuit of Eqs. ( 1) and (3), with Eq. ( 11) being the electric counterpart of (3). [ Equations (5)](#m5){ref-type="disp-formula"} and (9), conversely, represent the corresponding piecewise linear models.  Given a set of "genetic" parameters *α~i~*, *β~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 and quality factor *Q* = *ω*~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[@b35][@b36].  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 , with *χ* = (*m*, *p*) and: for the biological case, and *χ* = (*x*, *y*) with: for the electric case.  In absence of external sources \[i.e.,  for *S~m~*(*t*) = 0 or *V*(*t*) = 0\], Eqs. ( 12) -- (13) possess an equilibrium point, which we indicate by (, ) or (, ).  Such a fixed point is promptly calculated from 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 *γ*~0~ can be recast in the form , thus indicating that the divergence of *γ*~0~ for a vanishing denominator is approached when 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](#f2){ref-type="fig"} summarizes our results in the case of a fluctuating sinusoidal periodic source *V*(*t*) = *V*~0~sin(*ω~p~t*) + *η*(*t*), being *η* a gaussian noise term with zero average 〈*η*〉 = 0 and time uncorrelated (white) variance 〈*η*(*t*)*η*(0)〉 = 2 T L *δ*(*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](#f2){ref-type="fig"}) and coherent ([Fig.  2d--f](#f2){ref-type="fig"}) contributions separately, as well as their combination ([Fig.  2g--i](#f2){ref-type="fig"}).  In order to simulate a realistic set of electric parameters, we considered *L* = 0.01 H, *C* = 310^−13^ F, *R* = 2 Ω, *V~f~* = −0.6 V and \|*V~T~*\| = 600 V. This choice gives a resonant frequency *ω*~0~ = 0.58 MHz and a quality factor *Q* = 2887.  In the example reported here, the source is characterized by *V*~0~ = 1.6 V, *ω~p~* = 1.4*ω*~0~, and *T* = 10^−7^ W.

[Figure 2a--c](#f2){ref-type="fig"} displays the system dynamics when only incoherent noise *η*(*t*) is present, with temperature *T* = 10^−7^ 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 *ω*~0~ and damping *ω*~0~/*Q*, and its main effect is to produces a "colored" noise ([Fig.  2b--c](#f2){ref-type="fig"}).  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(*ωt*) is considered ([Fig.  2d--f](#f2){ref-type="fig"}), 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 [](#f3){ref-type="fig"}[](#f4){ref-type="fig"}[Fig. ( 5)](#f5){ref-type="fig"} have been investigated in the past, both experimentally[@b37][@b38][@b39][@b40] and numerically[@b41][@b42].  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](#f2){ref-type="fig"}, the power density spectrum (here is the Fourier transform of *y*(*t*)) is particularly simple, it only shows a strong peak at the input frequency *ω* = *ω~p~* and a rectified component at *ω* = 0, which originates from the rectifying action of the diode.  The very interesting dynamics is observed in [Fig.  2g--i](#f2){ref-type="fig"}, when both coherent and incoherent sources are simultaneously present.  In this case, a strong harmonic component at a new frequency, *ω* = *ω~p~*/2, appears, which represents a subharmonic of the input frequency *ω~p~*, and period doubling is observed in the dynamics ([Fig.  2h--i](#f2){ref-type="fig"}).  Quite remarkably, the strong noise resonance at *ω* = *ω*~0~ ([Fig.  2c](#f2){ref-type="fig"}) disappeared from the dynamics.  By comparing the spectra in [Fig.  2i](#f2){ref-type="fig"} and [Fig.  2f](#f2){ref-type="fig"}, we observe that the current power densities *P*(*ω*) at *ω* = *ω~p~* (i.e., the input pump frequency) and *ω* = 0 (i.e., the rectified component) are approximatively the same.  The strong subharmonic peak at *ω~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](#f2){ref-type="fig"}).  To further investigate this process, we calculated the intensity of the peak at *ω* = *ω~p~*/2 for different input frequencies *ω~p~* and input amplitudes *V*~0~, as reported in [Fig.  3](#f3){ref-type="fig"}.  The subharmonic peak at *ω~p~*/2 has a resonant-like intensity as a function of the pump frequency, having its maximum when *ω~p~* ≈ 3*ω*~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 *ω~p~*-dependent threshold.  From [Fig.  3](#f3){ref-type="fig"}, we observe that the efficiency of the subharmonic generation process strongly increases in the region where *ω~p~* ≈ *ω*~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](#f4){ref-type="fig"} summarizes our results by considering square pulses of time duration *δt*/*T* = 1/6, with frequency *ω~p~* = 2*π*/*T* and amplitude *V*~0~ = 2.2 V. The noise temperature is set to *T* = 5 · 10^−7^ 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](#f4){ref-type="fig"}), where we report a larger bandwidth of the peak at *ω* = *ω~p~*/2 with respect to the sinusoidal case ([Fig.  2i](#f2){ref-type="fig"}).  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~* = *y*(*nπ*/*ω~p~*), which are obtained by sampling the current at half of the input period *π*/*ω~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](#f4){ref-type="fig"}.  The figure shows that the "stochastic" period doubling observed in [Fig.  2](#f2){ref-type="fig"}--[4](#f4){ref-type="fig"} 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](#f4){ref-type="fig"}, 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](#f2){ref-type="fig"}--[4](#f4){ref-type="fig"}, we consider the piecewise linear model (9) and refer to [Fig.  5](#f5){ref-type="fig"}, 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](#f5){ref-type="fig"}), 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 ON→OFF shift) and **B** (OFF→ON 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](#f5){ref-type="fig"}).  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 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*: having introduced the two auxiliary functions *f*(*x*) = (*x* − *x~o~*)*θ*(*x* − *x~o~*) + *x~o~* and *h*(*x*) = 2*θ*(*x* − *x~o~*) − 1, with *x~o~* = *CV~f~* and *θ*(*x*) the Heavside step function.  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* = *η* = 0.  Due to the smallness of the subharmonic peak with respect to the input frequency at *ω~p~* ([Fig.  2i](#f2){ref-type="fig"}, [4c](#f4){ref-type="fig"}), we apply first order perturbation theory.  In particular, we indicate with \[*ξ*(*t*), *ζ*(*t*)\] the nonlinear, noise-free solution of Eqs. ( 15) for *η*(*t*) = 0.  Although we cannot write this solution in closed form, we know that this is strongly oscillating at *ω* = *ω~p~* in time, and lies in the curve **A-B-C-A** (or **B-D-E-B**) of [Fig.  5a](#f5){ref-type="fig"}.  In the presence of nonzero noise, we set *x*(*t*) = *ξ*(*t*) + Δ*~x~*(*t*) and *y*(*t*) = *ζ*(*t*) + Δ*~y~*(*t*).  By substituting the latter expressions into Eq. ( 15), and retaining only the first order terms in Δ*~x~*(*t*) and Δ*~y~*(*t*), we obtain the following dynamics: where: being *δ*(*x*) the Dirac-*δ* 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 *ξ*(*t*).  The first region is characterized by *ξ* + Δ*~x~* \> *x~o~*, which is represented by the segment **B**-**C**-**A** (or **B**-**D**-**E**) in [Fig.  5a](#f5){ref-type="fig"}, with points **B**, **A** (or **B**, **E**) and a small interval around them excluded.  In this case *θ* = 1 and *f* = *x*, which yields *h*′ = 0 and *f*′ = 1. [ Equations (16)](#m16){ref-type="disp-formula"} can be trivially solved in the frequency domain *ω* and the solution for the current (c.c.  stands for complex conjugate) is a damped harmonic oscillator, with maximum amplitude at *ω* = *ω*~0~ and bandwidth : being the Fourier transform of *η*(*t*).  The second region of the ON state, which is the most critical, is in the vicinity of point **A** and is represented by , with .  In this situation, the term Δ*~x~* leads the dynamics to continuously oscillate between the **A**-**B** (or **E**-**B**) segment (ON state) in [Fig.  5a](#f5){ref-type="fig"} 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 Δ*~x~* experiences a discontinuous dynamics originated from the term , with rapidly changing from zero to a small but finite value due to the noise fluctuations.  The current Δ*~y~*, conversely, evolves smoothly thanks to .  In this condition, equations of motion (16) become: The term *ζ*(*t*) appearing in the RHS of the first equation oscillates with frequency *ω~p~*, and this imposes a resonance condition with the oscillation frequency of Δ*~x~*(*t*).  In the Fourier domain, in fact, the only frequency admissible is the one where the term oscillates at the same frequency with Δ*~x~* and Δ*~y~*.  This condition imposes phase-matching *ω~p~* ± *ω* = *ω*, which is equivalent to *ω* = *ω~p~*/2.  The only oscillation that can be observed in the evolution of Δ*~x~*, Δ*~y~* is therefore at *ω* = *ω~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)](#m19){ref-type="disp-formula"} 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 *ω~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 current , we represent *δ*(*x*) = *H* · rect(*Hx*), being rect(*x*) = *θ*(*x* + 1 − 2) − *θ*(*x* − 1/2).  After taking the limit for *H* → ∞, we obtain the following solution: The spectral power of the current at *ω* = *ω~p~*/2 can be then obtained by combining in time the two expressions for appearing in Eqs. ( 20) and (18) for *ω* = *ω~p~*/2 with weights 1 − *α* and *α* 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 (*ξ*, *ζ*).  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*--factor , which here must be consider as a fitting parameter as it depends on the time spent by *ξ* and *η* in the different regions of the cycle.  After some straightforward algebra, we obtain the following expression for the current power density: being Ω = *ω~p~*/2*ω*~0~, \|*η~r~*/*R*\|^2^ the current noise spectral power at the peak frequency *ω~p~*/2. [ Equation (21)](#m21){ref-type="disp-formula"} has a nonlinear stochastic resonance at Ω ≈ 0.5, which predicts a strong noise amplification of the subharmonic *ω~p~*/2 when *ω~p~* ≈ *ω*~0~. In order to verify our theory, we calculated the behavior of *J~y~* versus Ω by a series of simulations with a sinusoidal input with a varying frequency *ω~p~*.  For each *ω~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 *α* = 0.51 and given by a nonlinear least square fit ([Fig.  6a](#f6){ref-type="fig"}).  We observed an excellent agreement between Eq. ( 21) and the results from numerical simulations, confirming the strongly resonant nature of the process. [ Equation (20)](#m20){ref-type="disp-formula"} allows also to interpret the quantitative differences between the power density spectra in [Fig.  2i](#f2){ref-type="fig"} and [Fig.  4c](#f4){ref-type="fig"}.  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 *ω~p~*, where each component resonate with its subharmonics and get amplified through parametric resonance.  This process is highlighted in [Fig.6b](#f6){ref-type="fig"}, where we superimpose different part of the current power density spectrum ([Fig.  6b](#f6){ref-type="fig"}) obtained from a numerical simulation with *V*~0~ = 1.6 V, *T* = 10^−6^ W and *ω~p~*/*ω*~0~ = 1.4.  As seen in the figure, the power density of the amplified noise near *ω~p~*/2 ([Fig.  6b](#f6){ref-type="fig"} solid green line) matches very well the spectrum of the input signal near *ω~p~* ([Fig.  6b](#f6){ref-type="fig"} solid red line).  In the case of short pulses, parametric resonance is triggered by a larger spectrum near *ω~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 *ω~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 *ω*~0~, is normally observed as the generation of different higher harmonics 2*ω*~0~, 3*ω*~0~, 4*ω*~0~, ...[@b43].  In this situation the generation of *ω*~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[@b44].  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 mid-infrared energy harvesting, where light-matter interactions are modeled by the same circuit of [Fig.  1b](#f1){ref-type="fig"}[@b45]. 

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. 

A. Fratalocchi acknowledges funding from KAUST (Award No.  CRG-1-2012- FRA-005).  For parallel computing, this research used the resources of the Primalight group at KAUST University. 

The authors declare no competing financial interests. 

**Author Contributions** A.F. and G.R. conceived the work and developed the theoretical analysis.  A.F. and G.R. reviewed and wrote the manuscript. 

![ Functional autorepressor scheme.\
(a) Genetic and (b) equivalent RLD circuit.  Panel (c) shows the nonlinear protein degradation function (solid line) and its piecewise linear approximation in the case of *χ~t~* = 0.2 (dashed line).](srep07088-f1){#f1}

![ Summary of numerical results for random sinusoidal sources.\
(a--c) Simulations with incoherent noise *V*(*t*) = *η*(*t*), (d--f) dynamics with coherent input *V*(*t*) = *V*~0~ sin *ω~p~t* and (g--i) results from the combination of coherent and incoherent signals *V*(*t*) = *V*~0~ sin *ω~p~t* + *η*(*t*).  Panels (a,d,g) show the time evolution of the input source *V*(*t*), panels (d,e,f) display the time dynamics of the current *y*(*t*) and panels (c,f,i) illustrate the current power density spectrum .  In the simulation we set *V*~0~ = 1.6 V, *T* = 10^−7^ W and *ω~p~*/*ω*~0~ = 1.4.](srep07088-f2){#f2}

![ Pseudocolor plot showing the behavior of the current spectral density *J~y~* at *ω* = *ω~p~*/2 for different input frequency *ω~p~* and various voltages *V*~0~.\
In the simulations we set the temperature to *T* = 10^−7^ W.](srep07088-f3){#f3}

![ Period doubling in the presence of short pulses with time duration *δt*/*T* = 1/6.\
(a--c) show the case for a fixed noise *T* = 5 · 10^−7^ W, while (d) plot the behavior of the minima of the current versus temperature *T*.](srep07088-f4){#f4}

![ System cycle in the (*x*, *y*) plane for a sinusoidal source (blue line) and square pulses of duration *δt*/*T* = 1/6 (green line).\
Panel (a) shows the case without noise, panel (b) displays the dynamics and with noise at *T* = 10^−7^ W. In the simulations we set *V*~0~ = 1.6 V. The arrows in Panel (a) indicate the direction of motion along the cycles **A** -- **B** -- **C** -- **A** and **B** -- **D** -- **E** -- **B**.](srep07088-f5){#f5}

![ Comparison between theory and numerical simulations.\
Panel (a) compares the normalized subharmonic current power density (circle markers) versus analytic theory (solid line) for different normalized frequencies Ω = *ω~p~*/2*ω*~0~. In the numerical simulations we choose *V*~0~ = 2 V, *T* = 10^−7^ W and change the input frequency in the range *ω~p~* ∈ \[*ω*~0~/2, 2*ω*~0~\].  For each *ω~p~* we averaged over 40 simulations to guarantee convergence of the mean value and standard deviation.  Panel (b) shows the behavior of the normalized power density spectrum in the surroundings of the input frequency *ω~p~* (solid red line) and near the subharmonic peak at *ω~p~*/2 (solid green line).  Spectra are calculated for a single numerical simulation with *V*~0~ = 2 V, *T* = 10^−7^ W and *ω~p~*/*ω*~0~ = 1.4.](srep07088-f6){#f6}
