Trade-offs and Noise Tolerance in Signal Detection by Genetic Circuits

Genetic circuits can implement elaborated tasks of amplitude or frequency signal detection. What type of constraints could circuits experience in the performance of these tasks, and how are they affected by molecular noise? Here, we consider a simple detection process–a signal acting on a two-component module–to analyze these issues. We show that the presence of a feedback interaction in the detection module imposes a trade-off on amplitude and frequency detection, whose intensity depends on feedback strength. A direct interaction between the signal and the output species, in a type of feed-forward loop architecture, greatly modifies these trade-offs. Indeed, we observe that coherent feed-forward loops can act simultaneously as good frequency and amplitude noise-tolerant detectors. Alternatively, incoherent feed-forward loop structures can work as high-pass filters improving high frequency detection, and reaching noise tolerance by means of noise filtering. Analysis of experimental data from several specific coherent and incoherent feed-forward loops shows that these properties can be realized in a natural context. Overall, our results emphasize the limits imposed by circuit structure on its characteristic stimulus response, the functional plasticity of coherent feed-forward loops, and the seemingly paradoxical advantage of improving signal detection with noisy circuit components.


Response of deterministic two-component networks to constant and oscillatory signals
For a general network of N components, our starting point is the kinetic equations for the evolution of each species dn i dt = J + i (n 1 , n 2 , . . . , n N ) − J − i (n 1 , n 2 , . . . , n N ), (1.1) where n i → macroscopic concentration of the ith molecular species, J + i → total synthesis rate (production flux), J − i → total degradation rate (relaxation flux), Following Paulsson [1,2] we define the elasticities, giving the relative change of the reaction fluxes for component i after a change in the j component. As in main text, we use bars to denote that concentrations and reaction fluxes are evaluated at equilibrium (t → ∞).
Another quantities of interest are the time scales τ i defined as [1,2] 1 SinceJ + i =J − i ≡J i , 1/τ i is usually taken as the degradation rate of the ith component assuming linear degradation. We note, however, that this very same definition can be applied to cases of non-linear degradation of a chemical species by a different one, when considering a quasi-steady state approximation as valid [3]. For instance, in the case of hyperbolic consumption of a molecular component n i by another species n j in an enzymatic reaction, the relaxation flux can be written asJ and one simply defines then an effective decay rate of the form 1 To quantify the responses of the network we use the susceptibilities s i , defined as the relative change in the concentration of component i at equilibrium, after a change in the input signal, (1.5) Hornung and Barkai [4] showed that the susceptibilities can be expressed in terms of the elasticities H ij by differentiating Eq. (1.1) at equilibrium with respect ton I . If we consider two isolated interacting molecular species i, j, the susceptibility of species i due to a change in species j is simply (1.6) We used double index notation here to denote that the susceptibilities s ij give the local response of the network element i due to its direct interaction with component j, and we term this as pairwise susceptibilities. The susceptibilities s i defined in Eq. (1.5) give the global response of element i due to a change in input signal, and contain the effect of signal propagation through all connected network components. It is straightforward to show following [4] that susceptibilities s i in a detection network of N components, can be expressed in terms of the pairwise ones solving the algebraic equation If the input is a constant signal, the amplitude of the output response is proportional to the global susceptibility s O (Fig. S1). If the input is an oscillatory signal, the amplitude of the response will also depend on signal frequency. To see this, consider a single species n with birth/death kinetics whose production rate is periodically forced dn dt = α · (1 + a n sin (ωt)) − δ n n. (1.8) This is a linear differential equation with a time-dependent term which can be readily solved, for instance, by a Laplace transformation. The result is n(t) = α δ n − c 1 · e −δnt + αa n δ n sin (ωt) − ω cos (ωt) ω 2 + δ 2 n . (1.9) The second term is a transient, therefore at long times n oscillates around the equilibrium value α/δ n . The sine and cosine terms can be written in a compact form by using the relation sin (ωt + arg(z)) = sin (ωt) |z| Re(z) + cos (ωt) |z| Im(z), (1.10) from which we get that the oscillating part has the amplitude A(ω) = a n α ω 2 + δ 2 n .
(1.11) Therefore, the squared amplitude decays as a function of the forcing frequency with a bandwidth ω BW = δ n .
If an oscillatory signal is propagated through a network of N components, the output will also oscillate around the stationary mean with the frequency of the input, and an amplitude dependent on this frequency and the network characteristics (time scales, susceptibilities, etc.). Since we are interested in the amplitude of the response around the equilibrium value, we linearize the kinetic equations defining the relative deviations from equilibrium as and then the amplitude of the relative deviations is given by To first order approximation, the relative deviations evolve according to the set of linear equations ∆n dt = M · ∆n + q(t). (1.14) The matrix M is the normalized Jacobian of the kinetic Eqs. (1.1), The vector q(t) contains the explicit time-periodic terms, which can be generally written as q i (t) = a i sin(ωt + φ i ). (1.17) Samoilov et al. [5] show that, at long times, the solution of the system of equations (1.14) oscillates with an stationary amplitude given by Therefore the relative deviations from equilibrium oscillate with a squared amplitude For the specific case of an oscillatory input signal n I , with dynamics described by Eq. (1.8), acting on two-component network of sensor n S and output n O species, the normalized Jacobian matrix reads 21) and the time-periodic vector q(t) is simply since we are considering propagation of a single oscillatory signal with amplitude a I and frequency ω I , impinging on a network at equilibrium. Then (1.23) We also note that the stability of the equilibrium state of the two-component module is given by the determinant and the trace of its Jacobian matrix. The Jacobian determinant is in this case given by, 24) and the trace by If Jac > 0 and T r < 0, the network steady state is stable. The self-elasticities are H ii = 1 if there is no autoregulation of the i network components, H ii < 1 for positive autoregulation and H ii > 1 for negative autoregulation. Therefore for the network to be in a stable equilibrium state we should keep 0 < H ii < 1 for positive autoregulation and s SO s OS < 1 in the case of positive feedback.
To check the validity of the linear approximations involved in the response to constant or oscillatory signals, we numerically solved the kinetic equations for a linear genetic cascade (see "Models" section) and obtained the amplitude of the output response for a step input signal of the type 26) or for an oscillatory input signal of the form (1.8). Within the linear approximation, the relative change in the output following a step signal of amplitude a I is given by the ω I → 0 limit of Eqs. (2)(3) in main text. This limit depends only on signal amplitude and output susceptibility, Thus, for all types of networks and small input perturbations, the relative amplitude A of the output response varies linearly with output susceptibility s O . This is shown in Fig. S1A for the genetic cascade. At a I = 0.01 (black circles), the response is exactly linear in the whole susceptibility range, while at perturbations higher than 10% of the equilibrium value ( a I = 0.1, red circles), it deviates from linearity, but increases monotonically with susceptibility. For an oscillatory signal the output will oscillate at the same frequency with an amplitude given to linear order by Eqs. (2)(3) in main text. In Fig. S1B we show that this approximation reproduces quite accurately the frequency dependent amplitude of a genetic cascade when signal changes up to 50% the equilibrium value (a I = 0.5, blue symbols).

Frequency detection properties of a linear cascade
As a measure of the frequency detection performance in different modules, we use the range of proper frequency transmission given by the system bandwidth, ω BW .
We employ the standard definition [16] i.e., the bandwidth is the range of frequencies where the squared amplitude of the response is above one half of its maximum value. For systems behaving as lowpass filters, such as the linear cascade, A max is given by the amplitude at ω I = 0.
Therefore ω BW is given by the solution of the equation (1.29) A look at Eq. (4) in main text shows that, for fixed time scales, the bandwidth is independent of the output susceptibility s O . How does it depend on the time scales of the cascade components? First, one sees that when all degradation rates δ i are similar, the larger the degradation rate (faster time scales) the larger the bandwidth, and viceversa (compare red curve with δ I = δ S = δ O = 2 and black curve, δ I = δ S = δ O = 1, in Fig. S2A). Moreover, being the product of three low-pass filters, the bandwidth is limited by the slowest time scale of the system. This is illustrated in Fig. S2A (blue curve), with the response of a linear cascade where the module components, sensor and output species, are fast (δ S = δ O = 2) but the input signal has a lower degradation rate ( δ I = 1) limiting the range of transmitted frequencies. In this case, the response is exactly the same if the slowest time scale is either in the sensor or the output species. Another feature of linear cascades is that adding successive layers slows down output response [6]. A single component with periodically forced production rate, Eq. (1.8), has a bandwidth given by its degradation rate, Eq. (1.11) and dot-dashed line/shaded region in Fig. S2B. If this single oscillatory component is taken as an input acting on a second element, the output response is given by the solution of Eq. (1.20) with which yields (

Simultaneous determination of bandwidth and susceptibility for two-component networks
In order to compare amplitude/frequency detection features of two-component network motifs in a consistent way, we fix the sign and susceptibility of the input/sensor interaction, which we denote as s * SI (s * SI = 2 for Figures in main text). Each two-component module is then characterized by the sign and ranges of its remaining interactions, given as pairwise susceptibilities s ij , and by the time scales of each component, given as the inverse of the degradation rates δ i . Throughout the paper we take similar time scales for each node in the network, including the signal, and fix δ i = 1. We also allow pairwise susceptibilities to vary in the range s ij ∈ [smin, smax]. We take smin = 0 and smax = 5 in main text.
Next we find, for the simplest two-component networks analyzed here (positive and negative feedbacks, positive and negative autoregulations, and coherent/incoherent feed-forward loops) which is the relevant interaction parameter determining the bandwidth of the oscillatory response. Plotting the bandwidth as a function of FS we obtain the black curves in Fig. 2 in main text. Next, we determine the output susceptibility range compatible with each FS value. This is easily obtained from Eq. (1) in main text. For instance, for negative feedback the output susceptibility is The amplitude/frequency detection features for negative(positive) autoregulation as a function of ARS are similar to those for negative(positive) feedback and we plot them in Fig. S3: increasing ARS in negative autoregulation improves bandwidth but decreases amplitude detection, while the opposite behavior is observed for positive autoregulation.
For modules with feed-forward connections between input and output elements, s OI = 0, the bandwidth is determined by the relative contribution of the second term in the numerator of Eq. (7) in main text, acting as a high-pass filter. Thus we define the relative strength (RS≡ |s OI /s O |) as the relevant interaction parameter for frequency detection in feed-forward loops. It is clear that the highpass filter term dominates when s 2 OI ≫ s 2 O , which is only possible for incoherent FFLs, where sign(s SI · s OS ) = sign(s OI ). This regime of high-pass filtering of an oscillatory signal effectively takes place for RS> 3.3 ( Fig. 3 in main text), meaning that low frequency oscillations are transmitted with an amplitude smaller than A max / √ 2 (strictly speaking, circuits in this regime are band-pass filters, since there is always a cutoff at very high frequencies due to the intrinsic time scales of the components, but we use the denomination high-pass in order to stress the difference with the low-pass regime at fixed time scales). The high-pass behavior of the incoherent FFL (illustrated with trajectories in Fig. 3F in main text) is shown in Fig. S4 for the whole range of input frequencies (grey symbols) at RS=6. In order to compare with the negative feedback case at the same susceptibility and similar frequency detection regime (FS=6), we also plot the negative feedback oscillatory response with black circles.
The relevant interaction parameters determining frequency detection, together with the corresponding susceptibility ranges for each of these parameter values, are provided in Table S1.

Amplitude of fluctuations within the linear noise approximation
For small number of molecules, as it is usually the case in signaling reactions inside cells, the proper mathematical framework to describe the network dynamics is the master equation for the evolution of probabilities of the number of molecular species [7]. The kinetic equations (1.1) describe only the evolution of average number of molecules [12] (related to the macroscopic concentrations n i by the "system size" or volume factor V as n i = V n i ) . The amplitude of the fluctua-tions around this average numbers at equilibrium, assuming small fluctuations, is usually obtained following two alternative routes 1) Starting from the master equation, there can be derived exact equations for the time evolution of first and second moments of molecule numbers [8,9]. Linearizing around the equilibrium state, one obtains a Lyapunov algebraic equation for the covariance of the fluctuations around steady state [9,10]. The linear approximation should be valid in principle whenever the equilibrium point is far from a bifurcation, and the size of the fluctuation is not too large. Numerical simulations show that for steady states far from bifurcations, this is a good approximation even for large fluctuations. For an analysis of fluctuations using the Lyapunov equation near critical points see [11]. With the same notation used in Section 1 for reaction fluxes and elasticities, the Lyapunov equation can be expressed as [1,2] where M is the normalized Jacobian matrix, whose elements are defined as in Eq. (1.16); σ is the matrix of normalized covariances, σ ij = (n i − n j )(n j −n j ) /n inj (note that σ ii is the squared coefficient of variation for fluctuations of the i component), and D the diffusion matrix whose elements depend on the reaction fluxes, system size and stoichiometric coefficients [1,2]. To solve the algebraic equation (2.34) we assume that each molecular reaction affects a single species (which is not always the case, but it holds for the genetic networks studied here) and thus D ij = 0, ∀j = i.
For an analysis of cases where the diffusion matrix is not diagonal but can be transformed to diagonal and the Lyapunov equation solved, see Ref. [9]. Moreover, if each reaction adds or removes a single molecule, the diagonal diffusion elements, which give the individual noise strengths D i , take the simple form [1] where V is the system volume. In practice we take a volume factor for each molecular species (V i ) as an effective parameter to change the noise strengths. Although we use Eq. (2.35) for the systems investigated here, all analytical results are expressed in terms of general noise strengths D i . For instance, if species i is produced as a burst of b molecules condensed in a single reaction, while decaying linearly with a time constant τ i (a situation commonly considered for mRNA transcription in prokaryotic gene expression [8]), φ its noise strength takes the form 2) A second derivation starts again from the master equation, and performs a system size expansion [12] which can be used to derive Langevin equations for the dynamical variables. The approximations involved in this procedure are discussed in Ref. [7]. Linearizing around the equilibrium states, one gets again a system of linear Langevin equations for the time evolution of the fluctuations. These define a multivariate Ornstein-Uhlenbeck process which can be solved by standard methods (Gardiner (1985) [13]), yielding again the Lyapunov equation (2.34).
Both approaches are equivalent and thus no additional approximations are involved. The advantage of using the Langevin formalism is that, as we see below, the power spectrum of the fluctuations can be readily obtained from the same matrix definitions as above.

Solution of the Lyapunov equation for two-component modules
It is instructive to consider first the solution of the Lyapunov equation (2.34) for a general system of two species, n 1 and n 2 . The Jacobian and diffusion matrices read (2.38) Species 1 may stand for mRNA and 2 for protein, for instance, and then we can study different forms of autoregulation. If they denote different proteins, one can treat all kind of feedback interactions between both species, whenever the condensed two-variable description is a proper approximation. The aim is to express the noise covariances in terms of the same deterministic quantities (susceptibilities, time scales and elasticities) used for the description of the amplitude and frequency response in Section 1. The result for the autocovariances is are the susceptibilities of species 1 and 2 respectively (in this case global and pairwise susceptibilities are the same), and the intrinsic noise prefactors are given by As discussed by Paulsson [1,2] T av can be interpreted as the total rate for returning to equilibrium (time scale for relaxation) after intrinsic fluctuations. For the case H 12 = 0, and with the noise strengths D i given by Eq. Next we solve the case studied here, that of a signal acting on a two-component module, with the Jacobian matrix given by Eq. (1.21). The only restriction is that signal dynamics is not affected by any of the module components. For a discussion on the effect of correlations among signal and downstream components see Ref. [14]. Proceeding as previously, solving the Lyapunov equation (2.34) and rearranging the results in terms of susceptibilities and time scales we arrive at the following expression for the autocovariance of the output fluctuations: where Jac is the two-component Jacobian determinant given by Eq. (1.21), T tot av the total relaxation time for intrinsic fluctuations, = 0 (only signal noise) has been analyzed by Hornung and Barkai [4], where they have shown that for the simple circuits studied here, positive feedback (or autoregulation) is the only architecture giving noise reduction with respect to the linear cascade at similar susceptibility values. The importance of taking into account all possible noise contributions is discussed in the next section, taking the coherent FFL as an example.
In Fig. S5A, the theoretical prediction for the output noise variance Eq. (2.44), expressed as a coefficient of variation, is numerically validated with Gillespie simulations [15] for three different networks: a linear cascade (red), a coherent FFL (blue) and an incoherent FFL (green). The noise strength of sensor and output species is fixed through their "system size" parameters V S = V O = 1, and signal noise is varied monotonically increasing V I . We note that for the same susceptibility the coherent FFL (blue) can be tuned to decrease total output noise with respect to the linear cascade.

Analysis of noise/susceptibility trading in coherent FFLs
From the results of Fig. 5B in main text and the previous Fig. S5A we conclude that coherent FFLs usually have better SNR in amplitude than linear cascades. In order to see how this is possible, we compare the total output noise for both networks at the same total susceptibility, s lc O = s cffl O ≡ s O (recall that input/sensor interaction s SI is also fixed). To simplify the analysis, we assume equal time scales and noise strengths in all components: Then substituting in Eq. (2.44) the output noise variance for the linear cascade can be written as where we have used s OS = (s O − s OI )/s SI for the sensor/output susceptibility in coherent FFLs. Then we see that the extra direct interaction s OI adds noise to the term accounting for signal noise propagation(third term in brackets), but reduces noise propagated through the sensor species (second term). We can calculate the conditions for which the effective noise reduction in the sensor species is larger than signal noise amplification. Expressing everything in terms of the relative strength (RS) defined above for feed-forward loops, we find that (2.48) for total noise reduction, which depends only on the direct input/sensor interaction, and not on total susceptibility. Thus, for s SI ≤ 2 we have noise levels in the coherent FFL lower than those of a linear cascade at equivalent signal sensitivity (recall that RS ≤ 1 for the coherent FFL, and that this analysis holds for equal noise strengths and time scales). This is illustrated in Fig. S5C, where we plot output coefficient of variation as a function of output susceptibility at s SI = 2 for the linear cascade (red) and the coherent FFL (blue). The rest of the interactions are varied sampling the interval s ij ∈ [0, 5] uniformly, as in Fig. 5 in main text.

Power Spectra and Correlations of Fluctuations
To analyze fluctuations in the frequency domain, it is useful to start from the Langevin approximation for the stochastic dynamics [7,12,17]. Linearizing around steady sate, one gets linear evolution equations for the fluctuations around equilibria. The dynamics for the normalized fluctuations η i = n i −n ī n i is given by where the diffusion matrix D is the same appearing in Eq. (2.34). Eqs. (3.49)-(3.50) define a multivariate Ornstein-Uhlenbeck process whose solution is well known [10,13]. The power spectrum matrix in stationary state, P η ij (ω) = η i (ω)η j (ω) follows as [13] P η (ω) = 1 2π We note that the correlations C ij = η i (t)η j (0) can be also obtained directly from the Langevin formulation using the regression theorem [13,17]  with σ being the variance matrix. The power spectra can also be calculated from the Laplace transform of the correlations as

Power spectra for two-component networks
Proceeding as for the noise variance, we solve first Eq. (3.51) for the simpler case of two interacting molecular species n 1 and n 2 (mRNA/protein, protein 1/protein 2, etc.) with Jacobian and diffusion matrices given by Eq. (2.38). After some algebra we rearrange terms as Similarly to the case of noise variance, the power spectrum for species n i can be decomposed into an intrinsic contribution whose magnitude at zero frequency is proportional to its noise strength D i , and an extrinsic contribution depending on the other species' fluctuations and susceptibility. Consider the case that n 1 acts directly upon n 2 with no feedback (H 12 = 0) (3.60) Then one sees that n 1 (input) has the spectrum of a pure Poisson or birth/death process, which is a low-pass filter with bandwidth ω BW = H 11 /τ 1 , whereas the n 2 spectrum (output) has an intrinsic part given by its own Poissonian filter with bandwidth H 22 /τ 2 , plus an extrinsic contribution which is the product of the input and output filters scaled by the static susceptibility [14,18]. Solving Eq. (3.51) for our input/sensor/output network with Jacobian matrix Eq. (1.21) is equally straightforward and gives Eq. (11) in main text for the output component. We recall that this gives the spectrum of stationary fluctuations after the output element reached steady state (for instance, following a change in signal amplitude). For an oscillatory signal with noise, the output spectrum P (ω) can be decomposed into an oscillatory part and a noise background [19], where the background term P f luc. (ω) is the power spectrum of the output fluctuations around the mean, Eq. (11) in main text, and the oscillatory component is given by Here A(ω I ) is the amplitude of the output oscillations, given to linear approximation by Eq. (2) in main text. The signal-to-noise ratio (SNR) can be calculated as [19] SN R f req (ω I ) = 2 lim δω→0 . (3.63) To check the validity of the above expressions we show the numerical output power spectra of three different genetic circuits in Fig. S5B. The theoretical background spectra given by Eq. (11) in main text are plotted in red (linear cascade), blue (coherent FFL) and green (incoherent FFL). The height of the peaks at the oscillatory input frequency ω I = 1, marked with crosses, are proportional to the amplitudes A 2 (ω I ). This is shown in the inset of Fig. S5B, where we plot the theoretical A 2 (ω I ) given by Eq. (2) in main text for the linear cascade (red line) and incoherent FFL (green line) and the numerical peak heights of those modules at different input frequencies (black circles). We note that the numerical power spectra are divided by the bin size ∆ω (the frequency resolution) to be independent of the total sampling time, and thus are power spectral densities (power per unit frequency). While the background term P f luc. (ω) is a power spectral density, the oscillatory term in Eq. (3.62) is the product of a δ function (with units of inverse frequency) and the total power under it, therefore in order to compare with the numerical power spectra the theoretical oscillatory power must be divided by the resolution ∆ω [20].

Analysis of frequency SNR for different two-component modules
Similar to the analysis of noise/susceptibility features for the coherent FFL above, we simplify parameters not depending on circuit architecture as much as possible, and assume again equal time scales τ and noise strengths D for all network elements. Since we are comparing SNR of different circuits at input frequencies where frequency transmission is maximum [A(ω I ) = A max ], this usually happens at ω I ∼ 0 (except for the incoherent FFL which can behave as a high-pass filter for oscillatory signals). Substituting the theoretical expressions for P f luc. (ω I ) and A 2 (ω I ) in Eq. (3.63) and taking the limit ω I → 0 we have the following expression for the SNR where Jac is given by Eq. (1.24). Then, for linear cascades and FFLs with no autoregulation, Jac 2 = τ 4 and the SNR reads As in the analysis of noise amplitude above, we recall that s OS = s O /s SI for the linear cascade, but s OS = (s O − s OI )/s SI for the FFL. Therefore, for the same output susceptibility s O , giving the same oscillatory amplitude at ω I → 0, Eq. (1.27), the coherent FFL (s OI > 0) will have a frequency SNR larger than the linear cascade. The reason, as in the previous analysis, is again the reduction in the contribution to background power spectrum of the noise propagated through sensor component. This is shown in Fig. S5D where P f luc. (ω I ) is plotted versus A 2 (ω I ) at maximum transmission for the same sampling used in Fig. S5B.
For the feedback interactions the situation is different, since Jac 2 = τ 4 (1 − F S) 2 [Jac 2 = τ 4 (1 + F S) 2 ] for positive(negative) feedback. Substituting in Eq. (3.64) we see that positive feedback always decreases SN R f req with respect to a linear cascade with the same oscillatory amplitude, while negative feedback has the opposite effect. Thus the frequency SNR behaves differently than the amplitude SNR. While amplitude SNR improves in positive feedback due to time averaging [4], frequency SNR deteriorates due to loss of stability of the steady state (in the limit F S → 1 the steady state becomes unstable and SN R f req → 0). Another observation for feedback regulation is that s 2 O Jac 2 = s OS s SI which is the susceptibility of a linear cascade; therefore the maximum SN R f req at any feedback strength can not be larger than the maximum SN R f req of the linear cascade with the same interaction strengths (as it is shown in Fig. 4C-D in main text).

Regimes of frequency noise filtering
Two terms potentially giving high-pass filter behavior appear in the fluctuation power spectrum, Eq. (11) in main text. One of them is induced by the intrinsic output fluctuations [first term in Eq. (11)], and may be at work whenever the second term in the denominator ∆(ω) is different from zero. This happens if there is a feedback interaction between the output and the sensor species. Intuitively, one could expect that a negative feedback, which accelerates dynamics [21] will be able to accelerate also the intrinsic fluctuations. This has been demonstrated theoretically [22] and experimentally [23] by Simpson and coworkers. They showed that negative feedback in simple gene autoregulation or in a two-gene system produced a shift of the power spectral density to higher frequencies. This result can be also deduced from Eq. (11) by explicitly calculating the bandwidth of the module noise power spectrum. Expressing the denominator ∆(ω) in the alternative form ∆(ω) = Jac 2 + Γω 2 + ω 4 (4.66) where Jac is the Jacobian determinant defined in Eq. (1.24) and Γ is defined as the bandwidth frequency of the module fluctuations is given by ω mod.f luc.
In a negative feedback, the elasticities H SO and H OS have opposite signs. Therefore the Jacobian determimant (1.24) increases and the parameter Γ decreases, which from Eq. (4.68) it is seen to increase the frequency of the fluctuations.
To find the regimes where the time scale of fluctuations can be well separated in frequency detection, we calculate the overlap between the bandwidth of A 2 (ω I ) and the bandwidth of the fluctuation spectrum, P f luc (ω I ). If this overlap is smaller than the oscillatory bandwidth, there is a range of output frequency response free of fluctuations. To quantify this filtering capacity, we define the filter range (FR) by where ω osc BW is the oscillations bandwidth, ovlp is the overlap between oscillations and fluctuations bandwidth, ω osc max is the maximum oscillation frequency transmitted, and ω f luc max the maximum fluctuation frequency transmitted. Therefore, if F R = 0, there is no filter range, since ovlp = ω osc BW . For 0 < |F R| < 1, there is a range in oscillation bandwidth free of noise. If |F R| > 1 the system works as a perfect filter, completely separating the time scale of transmitted oscillations and fluctuations for the whole circuit bandwidth. The sign of F R indicates if the system filters fast fluctuations (F R < 0, since ω f luc max > ω osc max ) or slow fluctuations (F R > 0).
As suggested by theory, we found only two type of circuits which could filter out fluctuations by frequency content. In Fig. S6A, we show how the FR of negative feedback circuits changes as a function of FS. We see that for F S > 5, (F R < −1, shaded region in Fig. S6A) a negative feedback module can behave as a perfect filter completely separating fluctuations from transmitted oscillatory frequencies (since fluctuations are shifted to frequencies larger than oscillatory response bandwidth, Fig. 6A in main text). On the other hand, incoherent FFLs may act as partial noise filters (0 < F R < 1, Fig. S6B) when input oscillations are at high frequencies, since in this case the oscillatory response is large and fluctuations are always slower (Fig. 6B in main text).

Optimal detection by two-component modules requires FFL interactions
In the main text we have analyzed the simplest two-component detection motifs (feedbacks, autoregulations and FFLs). Mixed modules, made of a combination of two or more 'pure' motifs, may be able to improve simultaneous amplitude/frequency signal detection exploiting the particular advantages of each simple architecture. To explore this possibility, we have generated a random sample (10 5 elements) of all possible two-component modules with a uniform distribution of interaction strengths in the range s ij ∈ [0, 5] as in Figures 2-5 in main text. Since we also fixed the direct input-sensor interaction (s SI = 2, we have two possibilities for the sensor-output interaction (positive or negative) plus three possibilities (positive, negative or zero) for the remaining four interactions (feed-forward, feedback and two autoregulations), giving 2 × 3 4 = 162 detection circuits. We calculated output susceptibility, bandwidth, amplitude SNR and frequency SNR for each motif in the sample. We first selected those sampling elements giving both susceptibility and bandwidth larger than the linear cascade limit (s O > 10, ω BW > ω LC BW = 0.51) and identified the specific circuits. For each of those modules we calculated the frequency of appearance in the selected (s O , ω BW ) region, divided by the frequency of appearance in the whole (s O , ω BW ) space (Relative frequency, shown as histograms in Figure S7A). In this way, relative frequencies larger than unity indicate that the specific motif is not only capable of simultaneous amplitude and frequency optimal detection, but also that optimal detection is robust, since for random susceptibilities there is more probability of appearance in the optimal region than in the rest of the space. We found, Figure S7A, that most of the circuits in the selected regime are pure coherent FFLs. The only mixed modules which appear with a significative frequency are incoherent FFLs with a combination of feedbak or autoregulations of different signs, but in a much smaller proportion than coherent FFLs. We have used the same procedure to obtain the relative frequency of motifs with optimal signal detection with noise, selecting those circuits with both amplitude and frequency SNRs larger than the maximum value allowed in a linear cascade configuration ( Figure S7B). Again, motifs involving coherent FFLs are necessary for optimal detection, although mixed FFLs with positive feedback or autoregulation appear with similar relative frequencies than pure coherent FFLs. This is in contrast with the noiseless situation, Figure S7A.

General models and numerical simulations
The three component networks (input, sensor and output) numerically simulated in Figures 2 and 3 of the main text follow the same general system of differential equations (ω I = 0 for constant input signals) with linear degradation rates δ X . The production terms F S and F O are given in terms of Hill functions f act , f rep for activation and repression of the form where k X is the corresponding activation/repression threshold and n the Hill coefficient [24]. The expressions for sensor and output production terms in the two-component modules simulated are listed in Table 1. For positive feedback and coherent FFLs we use a function with OR logic for production of sensor and output species respectively. The reason for this choice is that the linear cascade architecture is recovered in the limit of large activation thresholds for the additional interactions (f act (O, k O → ∞), f act (I, k I → ∞) for positive feedback and coherent FFL respectively). These expressions hold for independent regulation of a single gene by two species (no physical direct or indirect interaction between the two transcription factors), but the results discussed here are the same considering other types of logic production gates or additional cooperativity parameters.
In the numerical simulations, we fix the degradation rates (usually to unity) and equilibrium concentrations for each species (typically in the range 10 2 − 10 3 ). Since we control for interactions with predefined susceptibilities, we recalculate, either analytically or numerically, the basal production rates α X and activation/repression thresholds k X to achieve the desired equilibrium concentrations and pairwise susceptibilities. Hill coefficients n are varied between 2 and 6 depending on the wanted output susceptibility s O (the susceptibility is related to the steepness of an effective activation function [1,2] and thus for s O = 5, for instance, we should use n ≥ 6).
Deterministic kinetic equations have been solved using Fortran codes and a fourth order Runge-Kutta algorithm. For the stochastic simulations, we employed Gillespie's algorithm [15]. Propensities are given by the production and degradation terms in Eqs. (6.70) (where n X now represents number of molecules of species X) with a suitable rescaling of basal rates and activation/repression thresholds by a 'cell volume' factor V X : α ′ X ≡ α X · V X and k ′ X = k X · V X . The noise strength in each molecular species can thus be independently varied from the other components by its own volume factor. In the numerical simulations we change only the input signal noise keeping V S = V O = 1 (Fig. S5A).
The power spectral density has been numerically obtained by fast Fourier transform of the autocorrelation function (using FFT Matlab subroutines). We run large enough ensembles of stochastic trajectories with different initial conditions, ∼ 10, 000, to ensure convergence of the autocorrelations.

Model analysis from experimental data
In order to assess if the new detection properties discussed for feed-forward loops can be found in natural or synthetic genetic circuits, we analyzed three experimental works where simple mathematical models of incoherent and coherent FFLs can be fitted to experiments.
Kaplan et al. [25] have recently measured in vivo and at high resolution the production rate of the galETK operon for a wide range of concentrations of its two natural inducers, cAMP and galactose. Expression of GalE by these inputs is mediated by an incoherent FFL of type I (see Figure 7A), where the transcription factor CRP (activated by cAMP) activates both the galE gene and an intermediate transcription factor GalS, which is a repressor of GalE (whose activity is dependent on D-galactose). A second constitutive repressor, GalR, was removed from the original system. The GalE-∆GalR production rate was fitted by Kaplan et al. to a Hill function dependent on cAMP and galactose concentrations [25]. We used this model to obtain the amplitude response ( Figure 7B) and calculate the signal detection properties shown in Figure 7C. The dynamics of the incoherent FFL is given by the following equations: where C is activated CRP, G S ≡ GalS and G E ≡ GalE. We assume that cAMP activates CRP in a linear fashion [25,26]. Parameters showing a good fit to experimental data [25] are k C = k Gs = k g = 5 mM, h1 = h2 = h3 = 1. For simplicity we assume that all protein degradation rates are given by cell growth and division, δ = log(2)/τ div = 1, which implies an average cell division time of ∼ 40 min, and α C = δ, so that active CRP at equilibrium is equal to cAMP concentration. With these equations and parameters the GalE production rate (assuming CRP and GalS are at equilibrium) reproduce well the experimental input function of the galE-∆GalR system, compare Figure 7B in main text and Figure  2 in Ref. [25]. The negative autoregulation of GalS has not been included in the model since it has little effect on the active repressor level and non-monotonicity of the input function [25].
To understand the peaked shape of the production rate as a function of cAMP levels (amplitude filter) we fix the galactose concentration and examine how global and pairwise susceptibilities change with cAMP, Figure S8. The susceptibilities for this simple model are easily calculated from Eq.(7.72) : where the input I ≡ CRP, sensor S ≡ GalS and output O ≡ GalE and species concentrations are at equilibrium. We see that the global susceptibility s O (red line in Figure S8B) is positive and reaches its maximal values at the lowest cAMP levels, due to the contribution of the direct CRP/GalE interaction (black dashed line). The negative GalS/GalE susceptibility (green line) is very small at low cAMP concentrations since the amount of repressor activated by CRP is below threshold. Note that in Eqs. (7.73), s OI > 0 but it has the typical form of a repressor Hill function, decreasing with CRP concentration due to saturation of the GalE promoter. The s OS negative susceptibility, on the other hand, changes with GalS as an activator Hill function, and is maximal at high repressor levels. Therefore the global susceptibility decreases with cAMP and, when the inhibitory interaction dominates, becomes negative. The cAMP level of zero susceptibility corresponds to the maximum of the GalE production (asterisk and blue line in Figure S8A) since from this point GalE production starts to decrease. This situation is independent of the specific models and interactions, and is expected for general incoherent FFLs showing an amplitude filter regime. As another instance, we analyze the model proposed by Basu et al. [27] which qualitatively reproduces the experimental response of a synthetic band-detect circuit to an external acyl-homoserine lactone (AHL) signal, Figure S9A. The inducer AHL activates the expression of lambda repressor (CI) and Lac repressor (LacI). LacI expression is inhibited by CI, effectively forming an incoherent FFL starting with AHL/LuxR (input) and ending at LacI (output). A GFP reporter under control of LacI repressor monitorizes the system response, Figure S9A. The dynamical equations provided by Basu et al. [27] are where A ≡ external AHL concentration, R ≡ LuxR/AHL complex, C ≡CI, L ≡LacI and G ≡GFP. Parameters are those fitting the experimental data of BD2 strain [27], Similar to the galactose system, the GFP response behaves as an amplitude filter for external AHL concentration, Figure S9B. Around to the maximum, where the susceptibility of the output element of the incoherent FFL (LacI) is close to zero, Figure S9C, the system filters low frequency oscillations and noise (inset of Figure S9B). The high-pass frequency filter for an oscillatory AHL signal is illustrated in Figure S9D, compare with Figure 3F in main text.
Different sugar utilization genes in E. coli are induced by cAMP through a coherent FFL architecture [26]. The inducer cAMP, produced in the cell upon glucose starvation, activates the master regulator CRP. This controls the expression of intermediate transcription factors(sensors) which are activated by specific sugars and, together with CRP, regulate the response of several genes (outputs) involved in the metabolism of the corresponding sugar. As in the previous example for the galactose system, the production rate of both sensor and output genes has been experimentally measured for coherent FFLs controling E. coli response to four different sugar systems: arabinose, maltose, rhamnose and fucose. Of these, one of the simplest to analyze is the maltose utilization module, Figure 7B in main text. Production of the sensor species, MalT transcription factor, is only dependent upon cAMP/CRP and not on maltotriose [26]. However, MalT binds promoters of the malEFG, malPQ and malK operons as an oligomer stabilized by the inducer sugar [28,29]. CRP and MalT both act cooperatively to activate MalE and MalK production [30]. With these facts in mind, we proposed a simple mathematical model of a coherent FFL based on Hill type regulation functions [24] to fit the experimental production rates of Kaplan et al. [26]. We specialized to the malEFG and malPQ operons, since MalK is known to directly interact with MalT in a negative feedback loop [28]. The equations for the CRP/MalT/MalE coherent FFL used are where s is the inducer sugar(maltotriose) C stands for active CRP, M T for active MalT transcription factor and M E for MalE product. We fitted independently the M al T and M al E experimental production rates to the mathematical form given by Eqs. (7.75) using a nonlinear optimization method with constraints (fmincon subroutine in Matlab 7.3) restricting the parameters to positive values. We took β = 65 mM and α E = 50 mM from the experimental values of MalT and MalE maximal production rates, and fitted the thresholds K 1 = 5.6, K 2 = 7.6, K 3 = 7.1 mM −1 and the Hill coefficients h 1 = 1.5, h 2 = 2.4, h 3 = 1.3 to the experimental data for MalT and MalE (root mean square errors 0.06 and 0.07 respectively, which are within the experimental 10% mean relative error [26]. We allowed for a small basal (α T = 4.4 mM) activation of MalT by CRP, and set the degradation rate by cell growth and division δ = 0.65 · log(2) h −1 as measured from the experimental average division rate [26]. The production rate of MalE for the model given by Eqs. (7.75) is compared to the experimental production rate obtained by Kaplan et al. [26] in Figure 1. The simple model proposed is of course not unique, and different mathematical functions can be fitted to the experimental data. We also considered models where some MalE activation is possible either by CRP or MalT alone, and nonlinear posttrancriptional activation of MalT by the inducer sugar s. The quality of the fittings as measured by the root mean squared error(rmse), however, was the same than for the simple model (0.06 < rmse < 0.08) and the results for amplitude and frequency SNRs qualitatively similar to those shown in Figure 7E,F in main text.