Multi-Target Regulation by Small RNAs Synchronizes Gene Expression Thresholds and May Enhance Ultrasensitive Behavior

Cells respond to external cues by precisely coordinating multiple molecular events. Co-regulation may be established by the so-called single-input module (SIM), where a common regulator controls multiple targets. Using mathematical modeling, we compared the ability of SIM architectures to precisely coordinate protein levels despite environmental fluctuations and uncertainties in parameter values. We find that post-transcriptional co-regulation as exemplified by bacterial small RNAs (sRNAs) is particularly robust: sRNA-mediated regulation establishes highly synchronous gene expression thresholds for all mRNA targets without a need for fine-tuning of kinetic parameters. Our analyses reveal that the non-catalytic nature of sRNA action is essential for robust gene expression synchronization, and that sRNA sequestration effects underlie coupling of multiple mRNA pools. This principle also operates in the temporal regime, implying that sRNAs could robustly coordinate the kinetics of mRNA induction as well. Moreover, we observe that multi-target regulation by a small RNA can strongly enhance ultrasensitivity in mRNA expression when compared to the single-target case. Our findings may explain why bacterial small RNAs frequently coordinate all-or-none responses to cellular stress.


Introduction
Since their discovery more than thirty years ago, it has become clear that small RNAs (sRNAs) play a crucial role in regulating gene expression. sRNAs downregulate gene expression posttranscriptionally by pairing with target mRNAs through base complementarity. Complex formation with small RNAs competitively inhibits mRNA translation and/or induces mRNA degradation (reviewed in [1,2]). Importantly, the interaction between the sRNA and its target is non-catalytic in nature, since sRNA molecules are typically degraded along with their target, instead of being re-used to regulate other targets [3]. Such regulation is distinct from other post-transcriptional regulators such as RNAbinding proteins (RBPs) and microRNAs (miRNAs) which, in most cases, pass through multiple rounds of mRNA complex formation [4]. On the other hand, RBPs and miRNAs can competitively inhibit translation, and thus resemble sRNA action [5,6]. Notably, sRNA-mediated regulation conceptually differs from transcriptional repression, besides simply regulating a later step in protein biosynthesis: transcriptional repressors are typically present in vast excess over individual binding sites in the genome; thus, unlike sRNAs, the repressor pools are not depleted by binding to specific targets. Compared to other modes of regulation, sRNAs may thus confer unique dynamical features to gene expression.
The quantitative aspects of sRNA regulation were analyzed by various mathematical modeling studies, most of which assumed a purely stoichiometric mode of sRNA action [7,8,9,10,11,12,13,14]. Model-based analyses revealed that sRNAs binding their targets with sufficiently high affinity can establish a threshold-linear gene expression response at steady state [9]: the stoichiometric nature of sRNA action ensures that mRNA translation is almost completely suppressed as long as the sRNA concentration exceeds that of the mRNA (sub-threshold regime). In contrast, gene expression increases linearly with increasing mRNA transcription as soon as the sRNA is less abundant than the mRNA species (linear regime). Recent work revealed that miRNAs can generate similar threshold-linear behavior at the single-cell level [15]. Regulation by sRNAs has a clear signature not only for steady state expression but also during dynamic responses. For example, the system may initially need to eliminate excess unbound sRNA when reaching the linear regime; thus the kinetic profile of gene expression is characterized by a sharp delay [7,16].
Mathematical modeling studies revealed that steady state and temporal thresholds require that sRNAs bind strongly to their targets. Accordingly, some sRNA-mRNA complexes were shown to be stable and can be assumed to form irreversibly, as expected for sufficiently long RNA duplexes [17,18,19]. However, sRNA species vary as to the extent they complement their targets, and even in case of extensive complementarity base pairing may only occur over a limited region, the so called ''kissing complex'' [20,21,22]. In order to maintain their regulatory effects, many sRNAs require the presence of a specific RNA chaperone protein, Hfq, which is thought to melt inhibitory RNA structures and may have a bridging function in mRNA binding [22]. In fact, the presence of Hfq increases local concentrations of sRNAs and mRNAs that drastically enhances complex formation. For example, a 50-fold increase in mRNA-association rate in the presence of Hfq has been measured between rpoS mRNA and DsrA sRNA [23] or between ompA and MicA [24]. Thus, in living cells high affinity complex formation may be ensured by additional factors beyond simple base pairing. Notably, Hfq has also been implicated in protecting the sRNAs from degradation until paired and in recruitment of the degradosome to sRNA-mRNA complexes [22], suggesting intensive regulatory interactions between sRNAs and RBPs.
Small RNAs regulate a broad range of targets, many of which function during bacterial stress responses, e.g., by regulating outer membrane proteins, iron metabolism, quorum sensing or carbohydrate metabolism [22]. By being expressed only under specific conditions, sRNAs contribute to adaptation to environmental changes. On the other hand, constitutively expressed sRNAs appear to limit the expression of proteins that become toxic when expressed in large amounts [21,22,25]. Many sRNAs have the potential for interacting with multiple mRNAs, which allows them to coordinate whole physiological responses [26]. One of the first sRNAs observed to interact with more than a single mRNA target was DsrA in E. coli effecting translation of global transcription factors [27]. Another well-established sRNA example for controlling multiple targets is RyhB, which downregulates at least 18 operons encoding iron-using proteins [28]. Very recently, 16 regulated mRNA targets were identified for Spot42, a small RNA that modulates catabolite repression in E. coli [29]. For GcvB and RybB in Salmonella it has been demonstrated that a conserved domain of the sRNA is responsible for recognizing multiple mRNA targets encoding ABC transport systems or outer membrane proteins, respectively [30,31]. A recent systematic target profiling and validation approach revealed RybB and MicA to act each as global mRNA repressors during envelope stress [32]. In summary, many sRNAs (e.g. RyhB, Spot42, RybB, MicA, CyaR, DsrA, GcvB, OmrAB, RNAIII) turned out to act as multiple-target regulators [26,29,32]. It is very likely that multitarget regulation by a small RNA is a recurrent design principle adding a new layer of complexity for bacterial gene expression. Multi-target regulation by a single sRNA was analyzed in previous mathematical modeling studies (reviewed in [33]): These papers showed that a single sRNA is able to efficiently downregulate many targets at once ( [13,34]), and that mRNAs crosstalk between each other by sequestering the common sRNA regulator [8]. Mitarai et al. proposed the concept of hierarchical prioritization, where RNAs respond sequentially to decreasing levels of sRNA expression, depending on their affinity for the common sRNA regulator [11]. In this paper, we extend these previous theoretical studies on sRNA-mediated multi-target regulation, and re-evaluate the parameter requirements of hierarchical prioritization by analytical and numerical studies.
Multi-target regulation systems (Fig. 1A, top) are recurrent motifs in biochemical networks, and are commonly referred to as ''Single-input modules'' (SIMs). Depending on the kinetic parameters, SIMs can show two types of behavior ( Fig. 1A): (i) if all targets show similar affinity for the regulator, they may respond synchronously once a critical regulator concentration is reached. Such behavior has been termed ''synexpression'', and is thought to require fine-tuning of target affinities [35]. (ii) In the other extreme, low-affinity targets may be relieved from inhibition at lower regulator concentrations than high-affinity targets. Exper-imental studies confirmed that such hierarchical prioritization occurs in transcriptional responses to morphogen gradients [36]. In this work, we compare different SIM architectures and parameter regimes, and show that some intrinsically favor synexpression while others tend to promote hierarchical prioritization. Since coordinated regulation of functionally related proteins is thought to optimize cellular responses [35,37,38,39,40], we focused on a SIM architecture that favors synexpression even without fine-tuning of target affinities: sRNAmediated multi-target regulation. Using analytical approximations, we comprehensively characterize the parameter space, and show that robust synexpression at steady state is a general feature of the system, unless target affinities are extremely different. Numerical studies show that shared sRNAs are able to synchronize the temporal dynamics of mRNA expression as well; moreover, sRNAs can establish strong feedback regulation in larger networks. Taken together, we extend previous studies on mRNA crosstalk by sRNA sequestration [8], and show that sequestration effects may be beneficial in robustly synchronizing mRNA expression levels under various regulation regimes. Our findings suggest that post-transcriptional control may be more efficient in coordinating gene expression when compared to transcriptional modes of regulation.

SIM architectures differ in their ability to robustly coordinate gene expression responses
Single-input modules (SIMs), where a single regulator controls multiple targets (Fig. 1A, top) are recurrent motifs in biochemical networks. SIMs may result in synexpression or prioritization of the regulator targets ( Fig. 1A and Introduction). To better understand the requirements for synexpression and hierarchical prioritization, we analyzed the steady state behavior of simple SIM architectures by calculating the response to increasing regulator concentrations. We focused on regulation mechanisms controlling gene expression at transcriptional and post-transcriptional levels. Consider the following minimal model of cooperative transcriptional coregulation, where a transcriptional repressor R T inhibits transcription of the mRNA species m 1 and m 2 (1) The kinetic parameters V max , K D , n and k deg correspond to the maximal transcription rate, the dissociation constant of promoter binding, the Hill coefficient and the mRNA degradation rate, respectively. The steady state solution [m i ] = V max,i /k deg,i / (K D,i n +R T n ) implies that the repressor concentration for halfmaximal inhibition equals K D,i . Thus, the threshold repressor concentrations eliciting half-maximal mRNA downregulation depend in a linear manner on the repressor affinity for the promoter. Moreover, the distance between the mRNA doseresponse curves scales with a9 = K D,2 /K D,1 (cf. Fig. 1B). Therefore, the transcriptional switch easily accommodates expression prioritization and sequential regulation, as even moderate affinity differences: (i) allow the two mRNAs to respond at significantly different regulator concentrations. (ii) establish an exclusive expression regime of the low-affinity mRNA at intermediate regulator concentrations (depending on the steepness of the doseresponse curves). On the other hand, precise synexpression of the two mRNAs in the transcriptional switch can only be established for very similar affinities K D,i , and thus requires extensive parameter fine-tuning. In the following, we will discuss SIM architectures whose mRNA expression thresholds diverge in a lessthan-linear manner with affinity differences. We will argue that this property promotes robust synexpression, independent of precise kinetic parameter values.
In Eq. 1, repressor pools are not depleted by binding to promoter sites; this assumption is likely to be fulfilled in transcriptional networks, where abundant transcriptional regulators (cf., [41]) target only two specific promoter copies per cell, and are buffered by non-specific binding to the genome. To understand the potential role of regulator depletion, we analyzed an extended cooperative inhibition model, where the repressor concentration is no longer assumed constant (see Protocol S1). The corresponding dose-response curves in Fig. 1C reveal that the distance between repression thresholds scales less than linearly with the affinity ratio a: for example, the threshold repressor concentrations eliciting half-maximal mRNA downregulation differ by ,3-fold if the affinities differ eightfold (a = 8). Parameter insensitivity is even more pronounced in the lower part of the dose-response. In such a system, synexpression would require less extensive parameter fine-tuning, making it a more plausible mechanism for precise co-regulation. We conclude that some SIM architectures intrinsically favor synexpression while others tend to promote hierarchical prioritization. In the Supporting Information, we show that the partial threshold insensitivity observed in Fig. 1C is directly linked to regulator depletion and sequestration effects (see Protocol S1). Regulator depletion effects play a major role in post-transcriptional gene expression regulation [8,33]; we therefore decided to comprehensively analyze the parameter (in)sensitivity in a minimal SIM, where a small RNA inhibits two . The x-axis is normalized by the regulator concentration where the case a9 = 1 reaches 10% of maximal mRNA expression. (C) An extended cooperative repression model with regulator depletion is less sensitive to parameter alterations. The steady state doseresponse of a mass-action model explicitly describing three sequential binding steps of the regulator to the two target promoters was calculated numerically (see Protocol S1 for details). This system is a generalization of the minimal cooperative repression model (Eq. 1). High affinity regulator binding to mRNAs was assumed to ensure that regulator depletion effects are significant. The x-axis is normalized by the regulator concentration where the case a9 = 1 reaches 10% of maximal mRNA expression to allow direct comparison with panel B. (D) Mathematical model for sRNA-mediated co-regulation. The mRNA species (m 1 and m 2 ) are controlled by reversible binding to the shared sRNA (s), giving rise to inhibitory complexes, c 1 and c 2 . The monomeric species m 1 , m 2 and s are constantly synthesized, and all molecular species are subject to degradation. Throughout this work we assume that the inhibitory complexes c 1 and c 2 are formed with high affinity. doi:10.1371/journal.pone.0042296.g001 mRNAs (Fig. 1D). We focused on a system where the sRNA is codegraded with its targets, because such a mechanism favors synchronization and coupling effects when compared to a catalytic mode of action (cf. Discussion). It will be shown below that the thresholds of the sRNA system are even more robust than the ones in Fig. 1C.

Modelling the regulation of multiple mRNAs by a shared sRNA
Post-transcriptional co-regulation of multiple mRNAs by a shared small RNA was analyzed using the model depicted in Fig. 1D. The translation of the mRNA species (m 1 and m 2 ) is competitively inhibited by the sRNA (s). Additionally, the effective turnover rate of the mRNA pools is controlled by the sRNA, since the degradation rates of the inhibitory complexes, c 1 and c 2 , may differ from those of the free mRNA species. By employing the law of mass-action, the dynamics of post-transcriptional regulation can be described by the following set of differential equations (2 This system comprises 12 kinetic parameters, i.e., three synthesis rates (v syn ), five degradation rate constants (k deg ), and four rate constants describing complex association and dissociation (k on , k off ). Simplifying assumptions were made to derive analytical expressions for the dose-response behavior and parameter sensitivity of the system: It was assumed that the system has reached steady state, thus (initially) neglecting the temporal aspects of the response. As the steady state of the full system (Eq. 2) cannot be solved analytically, the affinity of the mRNA-regulator complexes, c 1 and c 2 , was additionally assumed to be very high (K d,i = k off,i /k on,i R0). As outlined in the introduction, this seems to be a reasonable assumption for many sRNA circuits. In case of strong binding, all mRNA will be bound into the inhibitory complexes (c 1 and c 2 ) as long as the sRNA is present in excess, i.e., In the following, we will focus on the opposite regime (v syn,m1 +v syn,m2 .v syn,s ) to understand principles governing the accumulation of free, translationally active mRNA.
Analytical approximation for strong post-transcriptional co-regulation At steady state, transcription balances RNA degradation, giving rise to the following set of algebraic equations (4) v syn,s~kdeg,s : ½szk deg,c1 : ½c 1 zk deg,c2 : ½c 2 &k deg,c1 : ½c 1 zk deg,c2 : ½c 2 v syn,m1~kdeg,m1 : ½m 1 zk deg,c1 : ½c 1 v syn,m2~kdeg,m2 : ½m 2 zk deg,c2 : ½c 2 Degradation of free sRNA is neglected, since the whole sRNA pool will be bound completely into the complexes c 1 and c 2 as long as v syn,m1 +v syn,m2 .v syn,s . Taking into account the binding equilibria of the complexes , one can solve for the steady state concentration of the free mRNA species (5) Importantly, the steady state of the 12 parameter system in Eq. 2 is fully determined by five lumped parameters (a, b, S, v syn,m1 / k deg,m1 , v syn,m2 /k deg,m2 ). The mRNA synthesis ratio b relates the production terms of the two mRNAs (6a). : Large values of a imply that: (i) the sRNA affinity for m 2 is larger than that of m 1 (K d,2 eff ,K d,1 eff ) and/or (ii) sRNA-mediated destabilization of m 2 is large relative to m 1 (k deg,c2 /k deg,m2 . k deg,c1 /k deg,m1 ). Thus, a.1 indicates that sRNA-mediated inhibition of m 2 is stronger than that of m 1 . For equal inhibition of both mRNAs (a = 1) the steady state solution simplifies to (7 Thus, accumulation of free mRNAs can only occur if the sum of the mRNA transcription rates exceeds the sRNA transcription rate (v syn,m1 +v syn,m2 .v syn,s ; cf. Eq. 3). Using numerical simulations of the full system (Eq. 2), we confirmed that Eqs. 5 and 7 approximate the steady state of the system very well as long as complex formation is sufficiently strong. Equation 7 describes a threshold-linear response with a switch from no expression to linear expression at S = 1. Thus, for equal inhibition of both mRNAs (a = 1), each mRNA behaves similar to the case where a single mRNA is strongly inhibited in a thresholdlinear manner by a small RNA [7,9] or a miRNA [15].

sRNA-mediated co-regulation can synchronize gene expression thresholds despite different affinities for individual mRNAs
To analyze the parameter sensitivity of mRNA co-regulation, it is instructive to analyze the remaining fraction of active mRNA (9) f~½ m i v syn,mi k deg,mi The free mRNA concentration in the absence of sRNA equals [m i ] = v syn,mi /k deg,mi . Thus, the remaining fraction of active mRNA quantifies the percent inhibition of each mRNA by the sRNA. Plotting f as a function of the stimulus S (Fig. 2) yields doubly normalized dose-response curves fully characterizing mRNA regulation in terms of the constants a and b (Eq. 5). In biological terms, the stimulus S may correspond to alterations in the sRNA transcription rate (v syn,s ). For equal inhibition strength (a = 1), the mRNAs show perfect co-regulation, since the dose-response curves completely overlap regardless of the synthesis ratio b (Fig. 2 A, B, C; Eq. 7). Unequal inhibition strength (a?1) results in preferential suppression of the high-affinity mRNA. Nevertheless, the gene expression thresholds of the system respond in a sub-sensitive manner to changes in a, thus supporting that the system favors synexpression: The stimuli of 10% or 50% mRNA expression (dashed horizontal lines in Fig. 2  A, B, C) can be used to quantify threshold positions in a robust manner. Even if strong affinity differences are allowed (a#100) these threshold measures differ less than two-fold between the two mRNAs in most cases, and the ratio never exceeds a factor of 5 in Figs. 2 A, B, C. Thus, the thresholds of the sRNA circuit are even less parameter-sensitive than the thresholds of the cooperative switch with regulator depletion (Fig. 1C). Therefore, repression of the two mRNA typically occurs at very similar sRNA concentrations. Moreover, an exclusive expression regime for the lowaffinity mRNA ('sequential regulation') can only be established for very large affinity differences (a$100).
In the context of stress responses, it is most likely the all-or-none transition from no expression to a significant level that matters. In the following, we will therefore use the stimuli of 10% mRNA expression to systematically analyze alignment in this lower part of the dose-response. Synchronous switching is comprehensively quantified in Fig. 2D using the logarithmic difference between the 10% stimulus levels (S10) of both mRNAs (10). threshold ratio~log 2 S10(m 1 ) S10(m 2 ) Alterations in the parameters a and b generally induce subsensitive changes in the threshold ratio, confirming that synchronous switching is a robust property of the system. Threshold desynchronization is restricted to the upper-right region of the heatmap where the inhibition strength a and the synthesis ratio b are large (Fig. 2C and D). In this non-robust regime a highly abundant weak affinity mRNA and a low abundance high affinity mRNA coexist. Once the sRNA concentration is too low to inhibit both mRNAs (i.e., if S,1), the high affinity mRNA pulls the regulator away from the low affinity mRNA. Thus, the low affinity mRNA starts to be freed from the sRNA, while the high affinity mRNA is still subject to inhibition, giving rise to two distinct thresholds (dashed lines in Fig. 2C). It should, however, be noted that the thresholds still desynchronize in a sublinear manner even in this most sensitive regime. Besides transcriptional induction of the sRNA, the circuitry in Fig. 1 may also be controlled by alterations in the expression of one or both mRNAs, raising the question of whether synchronous switching is maintained under these conditions. The dose-response curves in Fig. 2 are based on the assumption that a and b are constant, and therefore continue to hold for co-linear regulation of both mRNA transcription rates (v syn,s = const.; v syn,m2 = b N v syn,m1 with b = const.). Thus, a constitutively expressed sRNA synchronizes gene expression thresholds arising from transcriptional coregulation of the two mRNAs. Moreover, parameter-independent coordinated switching is also maintained for selective transcriptional regulation of one mRNA (e.g., v syn,s = constant; v syn,m1 = constant; v syn,m2 ?const; Protocol S2). This is due to the fact that high-affinity sRNAs serve as push-pull devices that couple the expression of an unregulated target mRNA with that of a regulated target mRNA. For equal affinity of both mRNAs (a = 1; Eq. 7), the strength of this push-pull effect can be quantified by calculating the gain of the free m 1 concentration with respect to the synthesis rate (11) G~d ln(½m 1 ) d ln(v syn,m 2 )~v syn,m 2 ½m 1 : d½m 1 dv syn,m 2~1 1z1=b Thus, for the case where m 2 transcription is regulated, the expression of m 1 can respond infinitely strong if the system is close to the threshold (S<1). The effect is enhanced if m 2 is more abundant than m 1 (b.1). The former condition explains our observation that gene expression thresholds are coupled, while the upper parts of the dose-response curves are less well coordinated.
Taken together, we comprehensively characterized the parameter space to show a shared sRNA establishes synchronous gene expression thresholds irrespective of the mode of transcriptional regulation in the circuit. Numerical simulations confirm the synchronization effect, but also reveal that the effect is less pronounced if the low-affinity mRNA binds too weakly to the sRNA (non-stoichiometric mode of inhibition): in this case, large sRNA concentrations are required to fully repress the low-affinity mRNA, and, as a result, the threshold ratio would be higher than predicted by the analytical approximation.
sRNA-mediated co-regulation may convert thresholdlinear mRNA dose-response curves into a very steep allor-none switch Our simulations in Fig. 2B revealed that the low-affinity mRNA can exhibit an extremely steep dose-response curve for b = 0.1 and a$10. To understand this effect intuitively, we will analyze the limiting case of strong affinity differences (a&1) in the following. The steady state solution of the low-affinity mRNA (Eq. 5) becomes (12) 1). The remaining fraction of active mRNA species (Eq. 9) is shown as a function of the normalized sRNA synthesis rate (stimulus S; Eq. 8). The dose-response curves of the two mRNAs completely overlap if both mRNAs are inhibited with the same efficiency (a = 1; Eq. 6; black line), while they diverge for increasingly different inhibition strengths (a = 10; a = 100). According to Eq. 5, the system behavior is solely determined by the lumped parameters a and b. (B) and (C) Steady state doseresponse curves of mRNA expression (same as in panel A) for unequal synthesis rates of both mRNAs (b = 0.1 or b = 10; cf. Eq. 4). (D) Gene expression threshold synchronization occurs over a broad range of kinetic parameters. The log2 threshold ratio (Eq. 10) was calculated for varying mRNA inhibition strength (a) and synthesis ratios (b) using Eq. 5. Thresholds were calculated using the 10% stimulus levels (S10, cf. horizontal green line in panels A-C), and divided to calculate the threshold ratio (Eq. 10). Blue areas in the heatmap reflect synchronous switching of both mRNAs. doi:10.1371/journal.pone.0042296.g002 The behavior of the low-affinity species m 1 is governed by the terms T1 and T2. For low stimuli, the term T1 dominates, and m 1 will not be significantly affected by sRNA-mediated regulation, .i.e., (13) In molecular terms, this regime can be considered to represent sequestration of the sRNA by the high abundance/high-affinity species m 2 . As the stimulus increases, sequestration is less efficient and m 1 is subject to inhibition as well (term T2 dominates). Thus, sequestration of the sRNA by m 2 shifts the onset of m 1 downregulation to higher stimulus levels. Since complete inhibition of both mRNAs is fixed to S = 1, a small fold-change in the input S may be sufficient to switch m 1 from low to high levels, indicating a strong increase in ultrasensitivity. Under the assumption that a&1, we can estimate the onset of m 1 downregulation by setting T1 = T2<0, and solving for S. One obtains (14) Thus, the larger the expression level of the high affinity mRNA (i.e., the smaller b), the larger the shift in the onset and the more pronounced ultrasensitivity of m 1 . For very low b, S onset approaches the stimulus level where the sRNA concentration equals the sum of mRNA concentrations (S onset <1), indicating a perfect all-or-none switch. We have thus shown that the switching performance of a threshold-linear sRNA system can be strongly enhanced by the presence of a high affinity competitor, which may either represent another target mRNA or a binding decoy.
Equation 14 represents the extreme case of very large affinity differences (a&1). In Fig. 2A, we see that enhancement of ultrasensitivity is already observed for moderate values of a and b, although the effects are less pronounced in this intermediate parameter regime. For example, the estimated Hill coefficient as a measure of ultrasensitivity already approaches n H <6 for a ten-fold affinity difference of the two mRNAs (a = 10 and b = 0.1). For comparison, one estimates n H <2 in the single-target case.
Very strong affinity differences between the two mRNAs may imply that the limit of strong binding may no longer hold for the low-affinity mRNA. We therefore performed numerical simulations under the assumption that the low-affinity mRNA no longer shows a threshold-linear response in the single target case (i.e., does not follow a purely stoichiometric mode of inhibition). Importantly, the presence of a high-affinity competitor enhanced ultrasensitivity under these conditions as well (not shown).
To conclude, co-regulation by a shared sRNA gives rise to particularly interesting dynamic behavior if a high-affinity, highabundance mRNA and a low-affinity/low-abundance mRNA coexist. In this scenario, co-regulation induces robust synexpression and also enhances ultrasensitivity of the low-affinity mRNA (Fig. 2B).

Regulation by a shared sRNA synchronizes the temporal induction of mRNA expression
Our analyses in Fig. 1 indicated that regulator depletion and stoichiometric coupling effects favor synexpression of two mRNAs. Since stoichiometric depletion effects should also occur under presteady state conditions, we numerically analyzed whether sRNAs synchronize the temporal dynamics of gene expression as well. For simplicity, the analysis was restricted to a reduced model previously introduced by others [9,13]  Here, mRNA-sRNA complex formation is assumed to be irreversible owing to high affinity and/or rapid decay of the complex. Numerical analyses of temporal RNA expression were performed to confirm that sRNA-mediated regulation synchronizes the temporal dynamics of gene expression (Fig. 3). We assumed step-like alterations of mRNA or sRNA synthesis rates; these simulations reflect initiation or termination of physiological stress condition where gene expression responses are regulated by altering mRNA or sRNA transcription [2,25]. Upon step-like upregulation of both mRNA synthesis rates, while sRNA synthesis rate remained constant (at a non-zero value), the multi-target system behaved similarly to the single-target case [7,8]; sRNA-mediated regulation establishes a delay that equals the waiting time required to deplete the initial sRNA pool (Fig. 3A,  left). The duration of the delay phase was exactly the same for both mRNAs despite differences in sRNA affinity and/or expression levels, implying that the onset of mRNA expression is always synchronized. This has important ramifications for synchronizing the expression of mRNAs whose synthesis rates increase in a temporally sequential manner (cf. Fig. 3E and below).
When comparing systems with and without sRNA-mediated regulation, we observe that sRNA regulation can synchronize mRNA accumulation beyond the delay phase as well (Fig. 3A,  left). In the sRNA-less system (dashed lines), the accumulation of mRNA is solely determined by mRNA half-lives [42], implying that the short-lived mRNA accumulates faster. In the sRNA system, the difference in mRNA response times can be less than the difference between the half-lives (Fig. 3A, solid lines). This synchronization effect is due to opposite effects on the two mRNAs: Accumulation of the short-lived, high affinity mRNA is slowed down relative to the sRNA-less system (red solid line), since sRNA-mediated degradation efficiently suppresses mRNA accumulation at early (but not late) time points. The long-lived, low affinity mRNA responds more gradually to sRNA-mediated degradation enhancement; the net effect is an increased apparent mRNA turnover rate and thus a shorter response time (blue solid line). We conclude that sRNA-mediated regulation promotes synchronization; it will be shown below that the synchronization effect is particularly pronounced if short-lived mRNA has higher affinity for the sRNA (as assumed in Figs. 3A and B).
Upon coordinated downregulation of both mRNA synthesis rates, we observe a beneficial effect of sRNA-mediated regulation as well (Fig. 3A, right): In the sRNA system, both mRNAs are degraded faster than expected from their half-lives, and this agrees well with the previous observation that sRNA regulation speeds up mRNA downregulation in the single target case [7,8]. In the multitarget case, this acceleration is beneficial for synchronization and we observe an overlay of three effects: Firstly, enhanced degradation by the sRNA reduces the absolute difference between mRNA time courses, and thereby establishes synchronization. Secondly, once the short-lived, high-affinity sRNA declined to low levels, the sRNA is redistributed to the low-affinity, long-lived mRNA, thereby establishing a coupling effect. Thirdly, mRNA expression is synchronously and completely shut-off once the free sRNA regulator starts to accumulate; this sharpens the decay when compared to exponential kinetics of the sRNA-less system, and further promotes synchronization.
The opposite regulatory case of step-like sRNA synthesis with constant synthesis of both mRNAs is shown in Fig. 3B. We observe similar effects to the case where the mRNA synthesis is regulated (Fig. 3A): Upon a step-like downregulation of the sRNA synthesis, a common delay period is established for both mRNAs (Fig. 3B,  left). However, the following accumulation of both mRNAs is solely regulated by their half-life times; thus, sRNA regulation is not beneficial, simply because sRNA is not expressed at all in this scenario. Upon the step-like upregulation of the sRNA synthesis rates both mRNAs are rapidly and synchronously down-regulated (Fig. 3B, right); this behavior is reminiscent of the corresponding scenario in Fig. 3A (right), and arises from degradation enhancement and a sharp shut-off once the sRNA is present in excess.
Summarizing both modes of regulation ( Fig. 3A and B), our numerical analyses suggest that three simple rules govern gene expression dynamics of the sRNA circuit: (i) Stoichiometric coupling establishes the same waiting time for all mRNAs upon a shift from no to significant mRNA expression. (ii) The subsequent increase in mRNA species is aligned when the sRNA is further expressed. (iii) The transition from high to no mRNA expression can be efficiently synchronized due to a combination of enhanced degradation and strong stoichiometric inhibition.
To confirm that synchronous rise and shut-down of mRNAs (properties ii and iii) are robust features of the system, we Synchronous temporal switching occurs mostly independent of kinetic parameters. Response time alignment in the sRNA circuit was compared to the corresponding sRNA-less system using the synchronization factor (Eq. 16). Synchronization factors smaller than unity indicate that sRNA-mediated regulation enhances synchrony in mRNA expression. Synchronization was analyzed as a function mRNA degradation rates and the relative inhibition strength of mRNAs (a0 = k on,1 /k on,2 ). The heatmaps show simulation results for a regulatory scenario where a coordinate, step-like increase (C) or decrease (panel D) in both mRNA synthesis rates is accompanied by a step-like change in the sRNA synthesis rate in the opposite direction; this counter-regulation assumption was necessary to eliminate delay phases, and to obtain the same steady states in sRNA and sRNA-less systems. Qualitatively similar results are obtained for coordinate regulation of mRNA (but not sRNA) synthesis rates. See Eq. 15 for model equations and Protocol S3 for kinetic parameters. (E) Regulation by a shared sRNAs establishes a delay frame for mRNA induction. Multiple mRNA synthesis rates sequentially increase at different times in a step-like manner (top). Synchronous accumulation of mRNA species occurs upon depletion of the sRNA pool. Late mRNAs whose synthesis starts after the delay period respond immediately due to lack of sRNA buffer (dashed orange line). See Eq. 15 for model equations and Protocol S3 for kinetic parameters. doi:10.1371/journal.pone.0042296.g003 systematically analyzed the parameter space for these scenarios ( Fig. 3C and D). We compare sRNA and sRNA-less systems by calculating a synchronization factor (16) synchronization factor~abs t50 m reg Here, t50 equals the time point where the mRNA concentration has reached half of the difference between initial and final steady states value, and thus measures the response time. The abbreviations m i reg and m i unreg refer to the sRNA circuit and to the corresponding sRNA-less system, respectively. The synchronization factor thus quantifies synchrony as the absolute difference between the response times of m 1 and m 2 , and compares sRNA and sRNA-less systems by taking the ratio. Low values of the synchronization factor (,1) imply that sRNA-mediated regulation promotes synchrony of mRNA expression.
For mRNA upregulation kinetics beyond the delay phase, we observe that sRNA-mediated regulation synchronizes expression over a broad range of parameter values (Fig. 3C). As explained in the context of Fig. 3A, the synchronization effect is particularly pronounced if the short-lived mRNA has higher affinity for the sRNA (upper-left and lower-right quadrants in Fig. 3C). Regulation by sRNAs is even more beneficial in coordinating mRNA downregulation, since strong alignment of mRNA time courses is observed for almost all mRNA half-lives and affinity differences (Fig. 3D). Adverse effects of sRNA regulation are restricted to a narrow corridor, where mRNA half-life times are nearly equal and the sRNA-less case already shows pronounced synchrony.
For sRNAs to be effective in synchronizing the temporal dynamics of gene expression in a physiological setting, they should be able to coordinate mRNAs that are induced early and late during cellular stress responses. Numerical simulations were thus performed under the assumption that mRNA synthesis rates sequentially increase in a step-like manner at different times (Fig. 3E). Indeed, sRNA-mediated regulation synchronizes the induction of constitutively transcribed as well as early and late induced mRNAs by establishing a fixed delay time (blue lines and orange solid line in Fig. 3E). However, the synchronization capacity of the system is limited to the time point where the excess sRNA pool is fully degraded; a very late mRNA whose transcription starts afterwards shows no synchronization effect but its accumulation starts without delay (orange dashed line in Fig. 3E).
Taken together, we conclude that shared sRNAs improve temporal synchronization of gene expression by a combination of stoichiometric inhibition and degradation enhancement effects.

Discussion
Single-input modules (SIMs), where a common regulator controls multiple targets, may exert synexpression or hierarchical prioritization (Fig. 1A). Here, we show that SIM architectures differ strongly in their parameter sensitivity, implying that subtle regulatory differences intrinsically favor synexpression or prioritization. We show that post-transcriptional co-regulation of two mRNAs by a shared sRNA is particularly parameter-insensitive, and thus strongly promotes synexpression without a need for parameter fine-tuning. We numerically confirmed that similar conclusions hold true for the case where more than two mRNAs are strongly regulated by a shared sRNA (not shown). Our results do not contradict previous reports showing that sRNAs can establish hierarchical prioritization [11]: however, we conclude that efficient prioritization would require extremely different affinities for individual mRNAs, at least in the case of high affinity binding.
Small RNAs were previously shown to optimize cellular stress responses by establishing threshold-linear behavior and temporal switching [7,9]. Cellular stress responses often involve the expression of multiple genes. Our present analyses suggest that sRNAs may be well suited to coordinate the expression of functionally related stress response genes: strong binding effectively suppresses all stress genes under normal conditions, while allowing for highly synchronized all-or-none induction above a critical threshold. Furthermore sRNAs are able to tightly synchronize the temporal induction of multiple genes. Coordinated expression is, however, most pronounced in the regime where the sRNA transcription rate is close to the sum of mRNA transcription rates (near-threshold regime). Thus coordinated expression over a broad expression range by a sRNA is a nonrobust property of the system and requires fine-tuning of parameter values. This implies that sRNA-mediated regulation may, for example, not be advantageous to coordinate the levels of constitutively expressed protein complex subunits, without finetuning of binding affinities to the sRNA. Previous modeling studies showed that sRNA systems may show large fluctuations in the near-threshold regime [10]. Assuming that noise mostly arises at the level of transcription and that sRNA-mediated regulation is fast, coupling will ensure that functionally related mRNAs fluctuate in a coordinated manner, thus preventing imbalances in cellular regulation.
Two central and experimentally testable predictions can be derived from our model: (i) mRNA pools targeted by the same sRNA may be coupled by sequestration effects, i.e., transcriptional upregulation of one target mRNA should affect the expression of other targets (Eq. 11). (ii) multiple target mRNAs should exhibit synchronous thresholds, independent of precise sRNA binding affinities. The first prediction is well supported by the existing literature, since mRNA targets of bacterial small RNAs and eukaryotic miRNAs appear to communicate with each other by sequestration effects [9,43]. The second model prediction concerning synchronized gene expression thresholds by posttranscriptional regulation could be verified experimentally by measuring the expression of small RNA targets with different affinities for the common small RNA regulator. Recently work by Hao et al [44] showed that affinities for small RNA regulators and therefore also repression strength can be gradually tuned by sequence alterations in the complementary region of the mRNA and the small RNA. The existence of synchronized temporal thresholds could be confirmed experimentally by measuring temporal kinetics of target expression responses, e.g., after induction of stress. Such temporal responses are already available for multiple targets of the small RNA Spot42 [29]. However, the proposed synchrony effect would require simultaneous analysis of all target mRNAs in the same biological sample.
Our analyses revealed that coupling of mRNA pools by sRNA sequestration effects is the molecular mechanism underlying robust synchronization of mRNA thresholds. Other post-transcriptional regulators like microRNAs or RNA-binding proteins are not degraded alongside with their targets, and thus employ a more catalytic mode of action when compared to bacterial sRNAs. Our unpublished observations suggest that coupling effects are also possible in catalytic system, although the phenomena are less pronounced than in Fig. 2. Accordingly, degradation of multiple proteins by a common degrading activity was recently shown to establish co-regulation of protein levels [45]. Similar to the concepts we present here, co-regulation was especially pronounced when the degrading enzyme was close to its maximum processing capacity (i.e., close to the threshold). We speculate that threshold synchronization by sequestration may synergize with other synchronization mechanisms. For example, it has been show that positive feedback amplifying an upstream regulator synchronizes downstream events in the cell cycle [46]. It is well possible that relatively weak positive feedback may be sufficient to perfectly align pre-synchronized sRNA multi-target systems.
In larger networks, sequestration-based coupling may profoundly affect the system dynamics, as previously discussed for signaling cascades [47,48,49,50,51,52]. We expect the same to be true for sRNA networks; for example, a hidden positive feedback and bistability can arise if our minimal reaction network is extended such that m 1 induces the transcription of m 2 ( Fig. 4A; green arrow) [48]: increasing m 2 levels reduce the amount of sRNA available for inhibition of m 1 , thus further enhancing m 2 expression in a positive feedback loop. Similarly, negative feedback can arise if m 1 represses the transcription of m 2 (red arrow). Strong negative feedback is known to render biochemical networks insensitive to perturbations [53,54,55]. Accordingly, our numerical simulations reveal that the negative feedback model exhibits adaptation which renders the steady state level of m 1 partially insensitive to the m 1 transcription rate (Fig. 4B). As expected, no adaptation is observed if sequestration-based feedback regulation by the shared sRNA is removed from the system (Figs. 4B, D, E). We conclude that the synchronization effects discussed in the context of Eq. 10 are sufficient to establish strong feedback regulation in larger posttranscriptional regulatory networks. Sequestration-based negative feedback may contribute to robustness observed in miRNA circuitries [56] which, in one case, was attributed to a single miRNA targeting multiple mRNAs [57]. Detailed simulations of complex post-transcriptional circuitries and quantitative experi-mental analyses are required to further define the roles of posttranscriptional feedback and feedforward loops in genetic robustness.

Supporting Information
Protocol S1 SIM architectures differ in their ability to robustly coordinate gene expression responses (related to Fig. 1 and incl. Fig. S1).

(DOCX)
Protocol S3 Regulation by a shared sRNA synchronizes the temporal induction of mRNA expression (related to Fig. 3).