Accurate encoding and decoding by single cells: amplitude versus frequency modulation
Gabriele Micali, Gerardo Aquino, David M. Richards, Robert G. Endres
aa r X i v : . [ q - b i o . S C ] M a y Accurate encoding and decoding by single cells:amplitude versus frequency modulation
Gabriele Micali , , , Gerardo Aquino , , David M. Richards , , and Robert G. Endres , , ∗ Department of Life Sciences, Imperial College, London, UK Centre for Integrative Systems Biology and Bioinformatics, Imperial College, London, UK Dipartimento di Fisica, Universit`a degli Studi di Milano, Milano, Italy ∗ E-mail: [email protected]
Abstract
Cells sense external concentrations and, via biochemical signaling, respond by regulating theexpression of target proteins. Both in signaling networks and gene regulation there are two mainmechanisms by which the concentration can be encoded internally: amplitude modulation (AM),where the absolute concentration of an internal signaling molecule encodes the stimulus, and frequencymodulation (FM), where the period between successive bursts represents the stimulus. Although bothmechanisms have been observed in biological systems, the question of when it is beneficial for cellsto use either AM or FM is largely unanswered. Here, we first consider a simple model for a singlereceptor (or ion channel), which can either signal continuously whenever a ligand is bound, or producea burst in signaling molecule upon receptor binding. We find that bursty signaling is more accuratethan continuous signaling only for sufficiently fast dynamics. This suggests that modulation basedon bursts may be more common in signaling networks than in gene regulation. We then extend ourmodel to multiple receptors, where continuous and bursty signaling are equivalent to AM and FMrespectively, finding that AM is always more accurate. This implies that the reason some cells useFM is related to factors other than accuracy, such as the ability to coordinate expression of multiplegenes or to implement threshold crossing mechanisms.
Author Summary
Signals, and hence information, can generally be transmitted either by amplitude (AM) or fre-quency (FM) modulation, as used, for example, in the transmission of radio waves since the 1930s.Both types of modulation are known to play a role in biology with AM conventionally associated withsignaling and gene expression, and FM used to reliably transmit electrical signals over large distancesbetween neurons. Surprisingly, FM was recently also observed in gene regulation, making their rolesless distinct than previously thought. Although the engineering advantages and disadvantages of AMand FM are well understood, the equivalent question in biological systems is still largely unsolved.Here, we propose a simple model of signaling by receptors (or ion channels) with subsequent generegulation, thus implementing both AM and FM in different types of biological pathways. We thencompare the accuracy in the production of target proteins. We find that FM can be more accuratethan AM only for a single receptor with fast signaling, whereas AM is more accurate in slow generegulation and with signaling by multiple receptors. Finally, we propose possible reasons that cellsuse FM despite the potential decrease in accuracy.
Introduction
Cells are exposed to changing environmental conditions and need to respond to external stimuli withhigh accuracy, e.g. to utilize nutrients and to avoid lethal stresses [1, 2]. To represent (encode) chemicalsin the environment, either ligand-bound receptors trigger chemical signals or ion channels allow entry ofsecondary messengers. These in turn activate transcription factors (TFs), which then regulate target-protein production (decoding). In eukaryotic cells, the conventional view is that the level of signalingwithin the cell directly encodes the external stimuli, with consequent gradual changes in the nuclearTF concentrations. This is effectively an amplitude modulation (AM) mechanism [3–10]. However,recent single-cell experiments also show pulsating signals [3, 11–14] and bursty entry of TFs into thenucleus [3, 10, 15–17], in close analogy to frequency modulation (FM). (Note that, although there isno modulation of an underlying carrier wave as in radio broadcasting [18], the AM/FM terminology iscommonly used in quantitative biology [10,15].) Although several hypotheses have been put forward, thebenefits and detrimental effects of either type of response remain largely unclear.There is experimental evidence that both types of modulation occur in gene regulation. For example,take the budding yeast
Saccharomyces cerevisiae . Under oxidative stress the nuclear concentration oftranscription factor Msn2 is proportional to the H O concentration, suggesting an AM mechanism (Fig.1A,B) [10]. However, in response to a calcium stimulus, Crz1, which is normally cytoplasmic, enters thenucleus in unsynchronized bursts, regulating at least a hundred target genes (Fig. 1C) [15]. The level ofstimulus affects only the frequency of bursts, not their amplitude and duration, which implies FM (Fig.1D and inset) [15,19]. Similarly, Msn2 and its homologue Msn4 exhibit FM under glucose limitation [10].Bursty FM is also found in bacteria and mammals, indicating that this is a general modulation schemeacross different cell types. For example, during energy-depletion stress, the bacterium Bacillus subtilis activates the alternative sigma factor σ B in discrete stochastic pulses, regulating around 150 downstreamgenes [20]. In addition, isoform NFAT4 in activated T-cells shows similar behavior [21].What are the relative benefits of AM and FM? One important issue is the susceptibility to noise, whichaffects the accuracy of sensing. For example, in broadcasting radio signals it is well known that FM is lessaffected by noise than AM. This is because noise mainly deteriorates the amplitude, which is where theinformation is stored in AM. A similar argument also favors action potentials in communicating neuronalsignals over long distances [22]. In contrast, it has been hypothesized that for other cell types, such asyeast, the bursty nature of FM may introduce more noise than AM, so that AM might be preferable(Fig. 2A,B) [15]. However, two recent articles (which we discuss below) disagree with this and suggestthat FM may still be more accurate [23, 24]. In addition, it is important to remember that there areoften other factors than noise minimization. For example, it has been suggested that, in situations wheremultiple genes need to be up or down regulated, FM can provide greater coordination and reliability (Fig.2C,D) [15, 19].Mora and Wingreen considered a model for a single receptor embedded in a cell membrane andcompared the noise in the output for two signaling mechanisms: continuous (CM) and bursty (BM)modulation [23]. In CM, the receptor signals continuously whenever a ligand is bound, whereas inBM the receptor signals for a short, fixed-sized burst only upon binding of a ligand. As we explainbelow, for multiple receptors these mechanisms become equivalent to AM and FM, respectively. Byconsidering integral feedback control, a common network for sensing concentration ramps and preciseadaptation [25–27], it was found that, for fast binding and unbinding, the noise in CM can be twice thatfrom BM, suggesting that FM leads to greater accuracy. Despite this unexpected result, there are twokey points that need further clarification. First, the response was only calculated to lowest order in thesmall-ramp parameter, thus neglecting any time dependence of the noise. Second, the derivation solelyrelied on the small-noise approximation, which might work well for fast signaling, but could be inadequatefor slow gene regulation.Similarly, Tostevin et al. found biologically relevant parameter regimes of promoter switching in generegulation in which an oscillating input can produce a more constant and hence less noisy protein outputlevel than a constant input with noise [24]. Although interesting, this model is restricted to decoding andlinear pathways, and requires fine-tuning. Its general applicability remains unclear, such as whether anoscillatory input signal can be replaced by random bursts and still remain more accurate than a constantinput. In fact, oscillating signals are well-known to maximize target responsiveness while bypassingdesensitization from constant signals [28]. They also globally entrain with its period robust to noise [29].Such oscillators are found in circadian clocks, segmentation clocks, cell cycle, p53 DNA repair pathways,as well as nuclear factor NF- κ B, epidermal growth factor ERK, cAMP and Ca signaling [17, 30–39].This leaves the question of the relative benefits of AM and FM (with respect to random bursts) largelyunanswered. Figure 1. Experimental evidence for amplitude and frequency modulation. (A and B)Example data showing amplitude modulation from [10]. (A) Single-cell nuclear localization of Msn2transcription factor in response to H O stress as a function of time. The stimulus profile (input) is astep change applied at t = 0 (inset) which applies to all figure panels. (B) Average time trace fordifferent concentrations of H O stress. (C and D) Example data showing frequency modulationfrom [15]. (C) Single-cell nuclear localization of Crz1 in response to calcium stress as a function of time,showing bursts of Crz1. (D) The average frequency of bursts against calcium concentration, showing anincreased frequency with increased concentration. (Inset) Burst duration distribution for low (bluebars) and high (red bars) concentration. Both histograms are well described by the Gamma distribution h ( t ) = te − t/τ b , with τ b = 70s (black solid line), demonstrating that pulse duration is independent ofcalcium concentration. Experimental data in arbitrary units (AU) of fluorescence.Here, we aim to investigate the advantages and disadvantages of CM and BM (AM and FM) forencoding and decoding of constant concentrations and ramps. To build intuition, we start with a singlereceptor/ion channel (CM and BM). We consider concentration sensing by a linear pathway, allowingus to gain exact results for different temporal regimes (as suitable for fast signaling and slow generegulation). To provide analytical results, we extend the single-receptor model for ramp sensing by Figure 2. Advantages and disadvantages of amplitude and frequency modulation.
AM maybe less noisy than FM (A,B), but FM may allow coordinated expression of many genes (C,D) [15, 19].(A) In AM, low/high stimuli result in low/high levels of transcription factor (TF) inside the nucleus.(B) In AM, different nuclear TF concentrations (blue and red curves) lead to gene expression ofproteins A and B (see orange and green promoter functions respectively) with variable ratios (order ofdot and square changes). (C) In FM, the stimulus strength only affects the frequency of bursts, nottheir amplitude. (Inset) Schematic of TF (purple dots) binding promoter P A of gene A (orange) andpromoter P B of gene B (green) with different binding strengths. (D) In FM, the nuclear TFconcentration is always the same during a burst, only the frequency of occurrence changes. As aconsequence, the protein ratio stays constant.Mora and Wingreen. First, we introduce an alternative mechanism to integral feedback, the incoherentfeedforward loop (another common pathway motif for ramp sensing and precise adaptation [40–42]).This allows us to generalize the model to more than one pathway. Second, by explicitly including thetime-dependence of signaling noise, we are able to provide first-order analytical results for the accuracyof ramp sensing. Taken together, a general principle emerges, favoring BM for fast signaling and CMfor slow gene regulation. Finally, we generalize to many receptors and ion channels, a far more realisticsituation for biological systems, allowing us to make connection with AM and FM. While we found thatAM is generally more accurate than FM, we speculate why cells may still utilize FM in certain cases ofgene regulation. Results
Cells sense external stimuli with cell-surface receptors and/or ion channels, which ultimately leadto changes in the concentration and dynamics of active transcription factors (TFs) inside the nucleus.Cells control the response at two different levels. Firstly, cell-surface receptors signal to regulate theactivity of TFs in the cytoplasm. Secondly, inportin and exportin regulate the entry of active TFs intothe nucleus, thereby regulating transcription (Fig. 3A). Here, we build a theoretical model that encodesinformation from an extra-cellular environment in an intra-cellular representation. We distinguish twoways of encoding this information: continuous modulation (CM) and bursty modulation (BM). Once theinformation is encoded, various proteins can act together to implement a response (decoding), involvingregulatory networks. To provide a general analysis for arbitrary noise we first address concentrationsensing in a simple linear pathway using the master equation. However, to derive analytical results forramp sensing and pathways with feedback we apply the small-noise approximation. We finally extendthese models to implement amplitude (AM) and frequency modulation (FM) for many receptors orion channels. Accuracy is assessed by comparing the protein output noise for the different modulationschemes, assuming that the signal is decoded by the average concentration.
Single-receptor/ion-channel model
Following Mora and Wingreen [23] we build a single-receptor model that implements CM and BM.We call the extra-cellular species c , which is encoded intra-cellularly by the signaling rate u . Assumingwe are in the fast diffusion regime in which each ligand molecule can bind the receptor only once, thereceptor can be in either of two uncorrelated states: on when bound and off when unbound. This allowsthe receptor activity, r ( t ), to be written mathematically as a binary response, which takes value 1 in the on state and 0 in the off state. The extra-cellular concentration c affects the unbound time intervals τ u ,such that the binding rate is given by h τ u i − = k + c ( t ), where k + is the binding rate constant. In contrast,the bound time intervals, τ b , are exponentially distributed random numbers with average h τ b i − = k − ,where k − is the unbinding rate constant, which is independent of the extra-cellular stimulus concentration(inset in Fig. 3A). As for ion channels, some are ligand-gated or regulated by receptors, while othersare voltage-gated and hence dependent on action potentials [43]. In all these cases the stimulus affectsthe opening or closing times. In CM downstream proteins are produced with a constant rate α duringeach on time interval, which leads to a signaling rate u CM = αr ( t ), while in BM ζ = αk − − molecules areproduced instantly at the moment of binding with rate u BM = ζ P δ ( t − t + i ), where t + i are the bindingtimes (Fig. 3B). This choice for ζ allows a meaningful comparison of CM and BM as both produce, onaverage, the same amount of intracellular species. General approach to concentration sensing exhibits two regimes of accuracy
In order to provide a general result for arbitrary input fluctuations, we write down the chemical masterequation. For simplicity, we only consider concentration sensing with c ( t ) = c , but the model can also beapplied to ramps. Furthermore, we assume a linear pathway in which the receptor/ion channel activity r directly regulates an output species with copy number n (with production rate u and degradation rate γ ) (Fig. 3C, left). Since the receptor/ion channel activity is a two-state system ( on / off ), there are tworesulting master equations for CM (one for each state) describing the probability of being in the on and NucleusVacuole Ca Crz1
Cell
Calcineurin Crz1p Ca
2+ inportinexportin
Figure 3. Schematic view of signaling and gene regulation. (A) Cartoon of
S. cerevisiae inpresence of extracellular calcium, considered a paradigm of bursty frequency modulation. Calciumenters through plasma-membrane ion channels and can be stored (released) in (from) vacuoles.Intracellular calcium activates calcineurin, which dephosphorylates Crz1p. Once dephosphorylated,Crz1 binds inporting Nmd5p and enters the nucleus. Exportin Msn5p subsequently removes Crz1 fromthe nucleus. Cytoplasmic calcium pulses may correspond to Crz1 bursts in the nucleus [15]. Red arrowsindicate movement while blue arrows stand for chemical signaling. (B) Single receptor/ion channelactivity, r ( t ) (blue line), depends on the concentration of extra-cellular stimulus c . The signaling rate u differs between continuous (CM) and bursty modulation (BM). In CM, u is constant rate α duringbound intervals, with p b the probability of being bound. In BM, ζ molecules are realized at the time ofbinding with τ bursts the duration between consecutive bursts (binding events). (C) Different regulatorynetworks. Linear pathway used for concentration sensing. Incoherent feedforward loop and integralfeedback control allow chemical ramps to be sensed. off states, i.e. p on ( n, t ) and p off ( n, t ): dp on ( n, t ) dt = γ ( n + 1) p on ( n + 1 , t ) + αp on ( n − , t ) + k + cp off ( n, t ) − ( γn + α + k − ) p on ( n, t ) , (1a) dp off ( n, t ) dt = γ ( n + 1) p off ( n + 1 , t ) + k − p on ( n, t ) − ( γn + k + c ) p off ( n, t ) . (1b)Note that α ≥ k − , so molecules are generally produced in the on state. In BM, instead, the masterequations which describe the probabilities p on ( n, t ) and p off ( n, t ) of having n proteins at time t , are givenrespectively by dp on ( n, t ) dt = γ ( n + 1) p on ( n + 1 , t ) + k + cp off ( n − ζ, t ) − ( γn + k − ) p on ( n, t ) , (2a) dp off ( n, t ) dt = γ ( n + 1) p off ( n + 1 , t ) + k − p on ( n, t ) − ( γn + k + c ) p off ( n, t ) , (2b)with burst size ζ a positive integer. We solve Eqs. (1a) and (1b) with generating functions and simulateEqs. (2a) and (2b) with the Gillespie algorithm (see Materials and Methods).Simulations via the Gillespie algorithm show different outcomes for fast (small-noise approximationlimit, Fig. 4A,B) and slow (Fig. 4C,D) dynamics of the receptor. For fast switching ( k + c , k − ≫ γ ), forboth CM and BM, the probability has an unimodal distribution (Fig. 4B). On the other hand, in theslow switching regime ( k + c, k − ≪ γ ), the probability distribution becomes bimodal for CM and unimodalwith a long tail for BM, leading to drastically increased noise (Fig. 4D). The unimodal distribution forBM, which is simply due to the use of infinitely short pulses, would become bimodal for finite widthpulses.In order to classify the different dynamics and to compare CM and BM for arbitrary noise, we requireinformation on the probability distribution of n output proteins. In particular, we study the average,variance and skewness (the latter is encoded in the third moment) of the distribution for both CM andBM. Constraining the average output of CM and BM to be the same (Fig. 5A,B), we identify tworegimes for fast dynamics: k + c < k − (Fig. 5C) and k + c > k − (Fig. 5D). Specifically, for k + c < k − ,BM is more accurate (Fig. 5C, inset), while CM is generally more accurate when k + c > k − (Fig. 5D,inset), except for minimal burst size ( ζ = 1). However, for slow dynamics (and hence large noise), CM isalways more accurate than BM. The study of the third moment shows that, for slow switching and hencebimodality, BM has large asymmetry (Fig. 5C,F).These observations can be explained as follows, using the fact that the receptor/ion channel canonly detect information from the extra-cellular environment during unbound ( off ) time intervals, as theextra-cellular stimulus only affects the binding rate (Fig. 3). For fast dynamics, the two regimes canbe understood by comparison with maximum-likelihood estimation (MLE), the most accurate strategyfor encoding [44]. MLE estimates the ligand concentration c ML = k − h τ u i − from the average unboundtime interval h τ u i . The bound time intervals are discarded as they only contribute noise [44]. BM, whichproduces fixed-size bursts at the times of binding, approaches MLE when the bound intervals are shorterthan the unbound intervals. In this case, the times of the bursts effectively estimate the unbound timeintervals (Fig. 3B, bottom) and BM is more accurate than CM. However, when the bound intervals arelonger than the unbound intervals, BM cannot estimate the unbound time intervals anymore and becomesless accurate than CM. Since CM produces protein during the bound intervals, it signals according tothe average receptor activity p b = h τ b i / ( h τ b i + h τ u i ) (Fig. 3B, top). Hence, CM effectively containsinformation on both bound and unbound intervals, and thus can still provide a reasonable estimate ofunbound time intervals. An interesting exception is αk − − = ζ = 1, for which BM becomes slightly moreaccurate than CM. In the latter case, since the rate of protein production during a bound interval in CMis very low, there is uncertainty as to whether CM actually produces protein or not, which reduces itsaccuracy. In contrast, for slow switching the burst size needs to increase since BM produces the samelevel of protein as CM. Hence, BM is always less accurate than CM, independent of whether bound orunbound time intervals are longer. While we analytically demonstrate the connection with MLE for fastdynamics in the next section, an extended discussion without comparison to MLE can be found in S1Text and S1-S3 Figs. Figure 4. The two regimes in the linear pathway model based on the master equation. (A-B) fast ( k + c = 20 s − , k − = 100 s − , γ = 0 . s − , α = 100 s − , ζ = 1) and (C-D) slow( k + c = 0 . s − , k − = 0 . s − , γ = 1 s − , α = 25 s − , ζ = 500) switching. (A,C) Protein number as afunction of time from Gillespie simulations for CM (blue lines) and BM (red lines). (B) The probabilitydistribution for n target proteins is unimodal for both AM (blue) and FM (red). (D) The probabilitydistribution is bimodal for AM (blue) and remains unimodal for BM (red) but with a long tail in theslow switching regime. Small-noise approximation to ramp sensing confirms two regimes for fast dy-namics
To further investigate fast dynamics, we extend an analytical model for ramp sensing in the small-noiseapproximation [23]. Considering the single-receptor described in Fig. 3A,B, we linearize the system byaveraging over a time much larger than the binding and unbinding times. We further assume exponentialdistributions for τ b and τ u so that (cid:10) ( δτ b ) (cid:11) = h τ b i = k − − and (cid:10) ( δτ u ) (cid:11) = h τ u i = ( k + c ( t )) − , where c ( t )increases only very slowly with time (see below). Hence, signaling noise arises in CM due to variablebound time intervals (ignoring stochastic production of protein during bound intervals), while in BM thebinding times (bursting times) vary. Without loss of generality, we set α = k − , which is equivalent to Figure 5. First three moments of the protein distribution in concentration sensing fromthe master equation.
Averages (A,B), variance (C,D), and skewness (E,F) as a function of thefrequency of binding events, f = k + c / (1 + k + c /k − ). (Insets) Magnification of small-noiseapproximation region (fast switching). Analytical results for CM (blue) and numerical results for BM(red) as function of the frequency of binding events (logarithmic scale). Two regimes are shown: k − = 10 k + c ( α = 100 s − , γ = 1 s − , ζ from 1000 to 1) (left column) and k − = 0 . k + c ( α = 10 s − , γ = 1 s − , ζ from 1000 to 1) (right column). Averages from CM and BM are constrained tobe equal, i.e. ζ = αk − − . Variances of CM and BM exhibit two different regimes for fast switching: for k + c < k − BM is more accurate than CM (inset in C), while for k + c > k − CM is generally moreaccurate (inset in D), except for ζ = 1. Third moments show that, for large noise, the probabilitydistributions become asymmetric. ζ = 1. Hence, as we show in S1 Text, for averaging time much longer than k − − and ( k + c ( t )) − , the0average and autocorrelation (variance) of u ( t ) are given by [23] h u ( t ) i = k + c ( t )1 + k + c ( t ) /k − , (3) h δu ( t ) δu ( t ′ ) i = g k + c ( t )(1 + k + c ( t ) /k − ) δ ( t − t ′ ) , (4)with h δu ( t ) i = 0 and g = ( (cid:10) ( δτ b ) (cid:11) / h τ b i = 2 CM1 + (cid:10) ( δτ b ) (cid:11) / (cid:10) ( δτ u ) (cid:11) = 1 + [ k + c ( t ) /k − ] BM . (5)Note that only the variance differs between CM and BM. In particular, in Eq. (5) the ratio k + c ( t ) /k − determines whether g is larger in BM or CM, which ultimately determines which scheme leads to theleast noise. BM has the lower noise only when k + c ( t ) < k − , i.e. when h τ b i < h τ u i . In particular, in thelimit of fast unbinding ( k + c ( t ) ≪ k − ), the signaling noise for CM is twice as large as for BM.Sensing temporal ramps, i.e. the change of concentration with time, is crucial for locating nutrientsand avoiding toxins. We start by considering a stimulus whose concentration is constant for t < t = 0: c ( t ) = ( c t < c + c t t ≥ , (6)for constants c and c with c t ≪ c . By applying Eq. (6) to Eqs. (3-5), the signaling rate can berewritten to first order as u ( t ) ≃ ( u + u t + δu t ≥ u + δu t < , (7)where u , u are functions of c and c , and δu is the noise described by h δu ( t ) δu ( t ′ ) i (given in S1 Text).The condition c t ≪ c is necessary so that u behaves linearly in time with u t ≪ u . Under thiscondition, the factor g BM of Eq. (5) becomes g BM ≃ k c k − | {z } g ∗ BM + 2 k c c k − t, (8)where g BM is given by g ∗ BM for a constant external concentration. We now assume that the extra-cellularstimulus is encoded in the signaling rate u which affects the production of two output proteins withconcentrations x and y . Specifically, we compare the output noise of x and y between CM and BM usingthe incoherent feedforward (Fig. 3C, middle) and integral feedback (Fig. 3C, right) loops. Incoherent feedforward loop
The incoherent feedforward loop is a network motif in which u directly affects two outputs x and y ,while y inhibits x (Fig. 3C, middle). The loop provides precise adaptation to a step-change in stimulusand can also be used for ramp sensing. Mathematically, we use the following two coupled stochasticdifferential equations, dxdt = k x (cid:18) f ( u ) g ( y ) − x (cid:19) , (9) dydt = u − k y y, (10)1where k x is the rate constant for production and degradation of x , while k y is the rate constant fordegradation of y , and f ( u ) and g ( y ) are specified functions. In order to have adaptation the variable y needs to evolve slower than x , which requires k x > k y . Here we choose f ( u ) = e bu and g ( y ) = e bk y y ,where constant b has units of time. This allows us to obtain an analytic solution (see S1 Text for details). Integral feedback loop
The integral feedback loop [23] is another network motif for precise adaptation and ramp sensing.Here, u affects x only (the main output), while x activates y and y inhibits x (Fig. 3C, right). Thegeneral equations for this model are given by dxdt = uf ( y ) − k x x, (11) dydt = k y ( x − , (12)where k x is the rate constant for degradation of x , k y is the rate constant for production and degradationof y satisfying k x > k y , and f ( y ) is a monotonically decreasing function of y . Specifically, we choose f ( y ) = e − by , where b is a dimensionless constant. This again produces an analytic solution (see S1 Textfor details). Small-noise approximation
To analytically solve Eqs. (9) and (10) for the incoherent feedforward loop, and Eqs. (11) and (12)for the integral feedback loop, we linearize these equations within the small-noise approximation, andassume that we are in the fast-switching regime. This allows us to find analytic solutions in a particulartime window and under certain conditions which we define in S1 Text. Specifically, for the incoherentfeedforward loop in the small-ramp regime, the average values of h u ( t ) i , h x ( t ) i and h y ( t ) i are determinedby the differential equations Eqs. (9,10). Although there are no steady states for ramps, h x ( t ) i and h y ( t ) i show time-dependent stable solutions h x ( t ) i = e bu ky , (13a) h y ( t ) i = u k y − u k y + u tk y . (13b)Introducing x = h x i + δx and y = h y i + δy into Eqs. (9) and (10) with subsequent linearization thevariance of the target-protein copy numbers can be derived (see Materials and Methods). To first orderin small-ramp parameters the variances of x for both types of modulation are D ( δx ( t )) E CM = ∆ (cid:20) g CM u − − c k + /k − k x + k y u + g CM (1 − k + c /k − ) u t (cid:21) , (14a) D ( δx ( t )) E BM = ∆ (cid:20) g ∗ BM u − − c k + /k − + 3 c k /k − k y u + (cid:0) − k + c /k − + 3 k c /k − (cid:1) u t (cid:21) , (14b)where ∆ = b k x e bu ky k x + k y )(1+ k + c /k − ) , and g CM and g ∗ BM are parameters discussed in Eqs. (5) and (8). Thecorresponding results for species y are provided in Eqs. (S56) and (S58), and plots for species x and y are shown in Fig. 6B,D.Consistent with the master equation, these results show again two regimes: ramp sensing is moreaccurate for BM if k + c < k − , while CM is more accurate otherwise. For a constant environment(zeroth-order with c = u = 0) the regime is largely determined by the factor g . If k + c < k − ,2 Figure 6. Two regimes in incoherent feedforward loop based on the small-noiseapproximation.
Output noise, i.e. relative variance of x (top) and y (bottom), as function of thenon-dimensional ramp time u t/u for k + c < k − i.e. h τ b i < h τ u i (left) and k + c > k − i.e. h τ b i > h τ u i (right). CM and BM are shown by blue and red lines respectively. (A,B) BM is more accurate than AMfor k + c = 10 s − and k − = 6 . × s − . (C,D) CM is more accurate then BM for k + c = 10 s − and k − = 6 . × s − . Remaining parameters: k + c = 10 s − , k x = 5 s − and k y = 10 s − . g BM = 1 + (cid:10) ( δτ b ) (cid:11) / (cid:10) δ ( τ u ) (cid:11) < g CM =1 + (cid:10) ( δτ b ) (cid:11) / h τ b i = 2 (Fig. 6A,B). This is because the variability of the bound intervals (cid:10) ( δτ b ) (cid:11) canbe eliminated in BM (but not in CM), and the unbound intervals are well approximated by the durationbetween bursts ( τ bursts in Fig. 3). For k + c ≪ k − , BM effectively implements MLE. In contrast, CMis more accurate for k + c > k − , where g CM = 2 and g BM > p b in Fig. 3B). These results also apply to ramp sensing since the accuracyof the downstream proteins (decoding) relates again to the factor g and hence to the ratio between thebound and unbound time intervals. The integral feedback loop in Eqs. (11) and (12) shows very similarbehavior (provided in S1 Text). The validity of our analytical results are confirmed by simulations of thestochastic differential equation for both pathways in S4 Fig. and S5 Fig.3 AM is more accurate than FM for multiple receptors/ion channels
To address the question of whether AM or FM is more accurate in encoding and decoding, we con-sider a straightforward generalization to multiple receptors (or ion channels) (see S1 Text and S6 Fig.for details). AM can be obtained by considering unsynchronized CM receptors. In contrast, the experi-mentally observed sporadic bursts of nuclear translocation [10, 23] and hence FM might be explained bysynchronized receptors that individually operate with BM.For N unsynchronized ( us ) receptors, the resulting average and variance of the signaling rate are h u ( t ) i usN = N h u ( t ) i and h δu ( t ) δu ( t ′ ) i usN = N h δu ( t ) δu ( t ′ ) i in terms of the single-receptor quantities.Consequently, the relative variance, given by the variance divided by the average-squared, scales with theinverse of the number of receptors ( N ). On the other hand, for N synchronized ( s ) receptors, the averageand variance of the signaling rate are given respectively by h u ( t ) i sN = N h u ( t ) i and h δu ( t ) δu ( t ′ ) i sN = N h δu ( t ) δu ( t ′ ) i . The relative variance is now independent of N . Hence, unsynchronized receptors (AM)have a reduction of noise by a factor N compared to synchronized receptors (FM).For slow dynamics, or fast dynamics with k + c > k − , CM is generally more accurate than BM (atleast for ζ > N receptors, AM is more accurate than FM by an even larger margin. Incontrast, for fast dynamics with k + c < k − , BM is more accurate than CM by at most a factor of 2 (Eq.(5)). But since AM is N times more accurate than CM, AM becomes more accurate for encoding thanFM for more than two receptors. Since our results from the previous sections show that larger signalingnoise leads to larger output noise, the same rule emerges for decoding.From a physical point of view, how can receptors act in a synchronized fashion? Receptors may becoupled by adaptor proteins or elastic membrane deformations, allowing them to act cooperatively [45,46].In conclusion, while for fast dynamics (small-noise approximation) BM can be more accurate than CMup to a factor of two, two receptors/ion channels are sufficient for AM to become more accurate thanFM. Since cells have thousands of receptors and ion channels, AM becomes the most accurate modulationscheme. Discussion
Cellular responses to extra-cellular stimuli involve both encoding the external stimuli by internal sig-nals (which is normally fast) and subsequently decoding via the regulation of protein levels (which isnormally much slower). The internal representation of the external signal falls into two broad categories:continuous/amplitude modulation (CM/AM), where bound receptors continually signal and the inter-nal concentration itself encodes the external signal, and bursty/frequency modulation (BM/FM), wherereceptors only signal when first bound and the signal is encoded in the frequency of peaks. Here, wecompared the output noise for both types of modulation in the presence of a constant and a linearlyincreasing (in time) external concentration. Besides considering a linear pathway, we compared two non-linear network motifs: the incoherent feedforward loop and the integral feedback loop. These loops areubiquitous in biological systems. For example, the incoherent feedforward loop is found in chemotacticadaptation of eukaryotes [40] and transcription networks in bacteria [41], and the integral feedback loopis found in chemotactic adaptation of bacteria [25, 47] and in eukaryotic olfactory and phototransductionpathways [27].We found that, for a single receptor or ion channel, BM can be more accurate than CM for fastdynamics. This situation can occur when the average duration of the active on state is shorter thanthe average duration of the inactive off state (Figs. 5 and 6). In this case, BM effectively implementsmaximum-likelihood estimation, the most accurate mechanism of sensing [44]. If instead more time is4spent in the on state, then CM is generally more accurate (except when the burst size is minimal, i.eone). The reason behind this effect, which we analytically prove within the small-noise approximation,is that CM has information about both the on and off states, whereas BM only knows when a switchfrom off to on occurs. As such, CM effectively implements Berg and Purcell’s classic result of estimatingligand concentration by time averaging [48] (see also Discussion in [44]). In addition, we found that forslow dynamics CM is always more accurate than BM, independent of whether more time is spent in the on or off states, due to increased burst sizes (Fig. 5). Taken together our results suggest that BM shouldbe more common in signaling pathways than in gene regulation.The generalization to multiple receptors/ion channels allows AM and FM to be compared. AM,which arises from unsynchronized CM receptors, has a reduced relative noise due to spatial averaging,while the relative noise in FM from synchronized BM receptors remains identical to the single-receptorresult. (Note the observed nuclear bursts of approximately constant amplitude and duration support ourFM mechanism [10, 15].) As a result, AM is always more accurate than FM for more than two receptors(S6 Fig.). Since cells have tens of thousands of receptors and ion channels, this implies that the reason thatFM is sometimes observed in real systems must have a different origin. At least three possibilities presentthemselves. Firstly, FM can help to coordinate gene expression [15,19], which is particularly useful whenhundreds of genes are controlled by a single transcription factor, such as during stress response [49–51].Secondly, FM can enhance co-localization of proteins inside the nucleus, providing another way to improvecoordination of gene expression [52]. Thirdly, as with oscillatory signals, bursts can be used to activatetranscription by threshold crossing [32] while avoiding desensitization [28]. This may then push the cellto differentiate into a new state (such as under starvation to initiate competence) [53, 54]. It is alsoworth noting that by using seemingly redundant isoforms (such as NFAT1 and NFAT4 during an immuneresponse), AM and FM can be combined to enhance temporal information processing [21].While providing intuitive insights, it is clear that our models are highly oversimplified versions ofsignaling and gene regulation in actual cells. One of the main reasons for this is that we used idealizeddelta-functions as pulses in BM (and hence in FM). However, for example, in the calcium stress-responsepathway in Saccharomyces cerevisiae (Fig. 3A) nuclear bursts of Crz1p are on average two minutes long(Fig. 1D, inset). Most likely cytoplasmic calcium spikes determine the nuclear bursts (Elowitz, personalcommunication), but since the mechanism of calcium spiking remains poorly understood, such burstsare difficult to model. A further limitation of our models is that bursts only relate to translocation,whereas additional bursts may occur further downstream during transcription [55] (e.g. due to promoterswitching [24]) and translation [56]. Future models may need to include these details.Our models suggest further experimental investigation in multiple areas. Firstly, the distribution ofburst duration affects factor g (Eq. (5)), so that g = 2 in equilibrium for a single-step process and poten-tially g < -independent transcription factors (as well as Ca -dependent co-regulated genes) arewarranted in order to verify coordination of multiple genes [15]. Finally, to see if bursts help jump startnew cellular programs ( i.e. transition into a new “attractor”), global changes in gene regulation can bemonitored.A general understanding of FM may help prevent developmental defects and human diseases. Indeed,several biomedically relevant transcription factors, such as NF- κ B, p53, NFAT and ERK, show oscillatory5pulsing or random bursting [16,17,33–36,54]. In fact, the destabilization of regulatory circuits can underliehuman diseases: studies suggest that the coordination of gene expression could be critical in maintainingthe proper functioning of key nodes in such circuits. For example, the NFATc circuit is cooperativelydestabilized by a 1.5-fold increase in the DSCR1 and DYRK1A genes, which reduce NFATc activityleading to characteristics of Down’s syndrome [16, 60]. However, ERK pulses are regulated by both AMand FM with the same dose dependence, and it remains unclear how they affect cell proliferation andthe relevance to cancer [36].Broadly speaking, temporal ordering (regularity or periodicity) serves at least two roles in livingsystems [61]: extraction of energy from the environment and handling of information. While the first roleis well studied in terms of molecular motors at the single-molecule level, the second role is intellectuallymore difficult to understand as it requires a broader, more global understanding of cells. We believethat future work that combines single-cell experiments with ideas of collective behavior and engineeringprinciples is most likely to be successful.
Materials and Methods
Master-equation model for concentration sensing
The master equations for continuous modulation (CM), Eqs. (1a) and (1b), can be solved at steadystate using generating functions. In particular, we derive the first three moments of the probabilitydistribution using the general model in [62]. When the system is in the on/off state, the production rateof species x is α on/off . The degradation rate γ is independent of the state of the system. The probabilitydistribution of n target proteins at time t is then described by dp s ( n, t ) dt = γ ( n + 1) p s ( n + 1 , t ) + α s p s ( n − , t ) + k ¯ s p ¯ s ( n, t ) − ( γn + α s + k s ) p s ( n, t ) , (15)where ¯ s = off (on) when s = on (off). By defining the generating functions G s ( z ) = ∞ X n =0 p s ( n ) z n , (16)and using Eq. (15), a solution for G s ( z ) can be found, which then readily gives the moments of p ( n, t ).In particular, the variance and skewness are given by (cid:10) δn (cid:11) = X s ( ∂ z z∂ z G s ( z )) (cid:12)(cid:12) z =1 − h n i , (17) (cid:10) n (cid:11) = X s [ ∂ z z∂ z z∂ z G s ( z )] (cid:12)(cid:12)(cid:12) z =1 . (18)Full details are given in S1 Text.In order to solve the master equation for bursty modulation (BM), Eqs. (2a) and (2b), we use theGillespie algorithm [63]. If the system is in the on state with n proteins at time t , it can either switch tothe off state with transition rate given by k − / ( k − + γn ) or else remain in the on state and lose a proteinby degradation. If instead the system is in the off state with n proteins at time t , it can either switchto the on state with switching rate k + c / ( k + c + γn ) and, via a burst, increase its number of proteinsto n + ζ , or again remain in the same state and loose a protein by degradation. The time step betweenreactions, δt , is chosen from an exponential probability distribution λe − λδt , with λ equal to the total ratethat at least one reaction occurs.6 ODE models for ramp sensing
The following method applies to both the incoherent feedforward and the integral feedback loop. Tosolve the ordinary differential equations (9-12) we linearize around stable solutions, x ( t ) = h x ( t ) i + δx and y ( t ) = h y ( t ) i + δy , and assume that small δu leads to small δx and δy . Note that when sensinga gradually changing ramp, h x ( t ) i and h y ( t ) i are not steady states. Defining X = [ x ( t ) y ( t )] T we canrewrite these equations as dX ( t ) dt + M X ( t ) = (cid:20) w δuz δu (cid:21) , (19)where the matrix M and the constants w and z are defined in S1 Text. Analytic solutions are onlyavailable when M is time-independent. As shown in S1 Text, Eq. (19) can be solved and written as anintegral, which can then be evaluated with, for example, Wolfram Mathematica 8. Acknowledgments
We thank Guido Tiana for helpful discussions, and Michael Elowitz for providing details on his data.GM and RGE acknowledge funding from ERC Starting-Grant N. 280492-PPHPI. GA was supported byLeverhulme-Trust grant N. RPG-181 and DMR was supported by BBSRC grant N. P34460.
References
1. Perkins TJ, Swain PS. Strategies for cellular decision-making. Mol Syst Biol. 2009 Nov;5:326.2. Kussell E, Leibler S. Phenotypic diversity, population growth, and information in fluctuatingenvironments. Science. 2005 Sep;309(5743):2075–2078.3. Behar M, Hoffmann A. Understanding the temporal codes of intra-cellular signals. Curr OpinGenet Dev. 2010 Dec;20(6):684–693.4. Black JW, Leff P. Operational models of pharmacological agonism. Proc R Soc Lond B Biol Sci.1983 Dec;220(1219):141–162.5. Mackeigan JP, Murphy LO, Dimitri CA, Blenis J. Graded mitogen-activated protein kinase activityprecedes switch-like c-Fos induction in mammalian cells. Mol Cell Biol. 2005 Jun;25(11):4676–4682.6. Yu RC, Pesce CG, Colman-Lerner A, Lok L, Pincus D, Serra E, et al. Negative feedback thatimproves information transmission in yeast signalling. Nature. 2008 Dec;456(7223):755–761.7. Bosisio D, Marazzi I, Agresti A, Shimizu N, Bianchi ME, Natoli G. A hyper-dynamic equilib-rium between promoter-bound and nucleoplasmic dimers controls NF- κ B-dependent gene activity.EMBO J. 2006 Feb;25(4):798–810.8. Werner SL, Kearns JD, Zadorozhnaya V, Lynch C, O’Dea E, Boldin MP, et al. Encoding NF- κ Btemporal control in response to TNF: distinct roles for the negative regulators I κ B α and A20.Genes Dev. 2008 Aug;22(15):2093–2101.9. Giorgetti L, Siggers T, Tiana G, Caprara G, Notarbartolo S, Corona T, et al. Noncooperativeinteractions between transcription factors and clustered DNA binding sites enable graded tran-scriptional responses to environmental inputs. Mol Cell. 2010 Feb;37(3):418–428.10. Hao N, O’Shea EK. Signal-dependent dynamics of transcription factor translocation controls geneexpression. Nat Struct Mol Biol. 2012 Jan;19(1):31–39.711. Paszek P, Jackson DA, White MR. Oscillatory control of signalling molecules. Curr Opin GenetDev. 2010 Dec;20(6):670–676.12. D’Andrea P, Codazzi F, Zacchetti D, Meldolesi J, Grohovaz F. Oscillations of cytosolic calcium inrat chromaffin cells: dual modulation in frequency and amplitude. Biochem Biophys Res Commun.1994 Dec;205(2):1264–1269.13. Berridge MJ. The AM and FM of calcium signalling. Nature. 1997 Apr;386:759–760.14. Berridge MJ, Galione A. Cytosolic calcium oscillators. FASEB J. 1988 Dec;2(15):3074–3082.15. Cai L, Dalal CK, Elowitz MB. Frequency-modulated nuclear localization bursts coordinate generegulation. Nature. 2008 Sep;455(7212):485–490.16. Crabtree GR, Graef IA. Bursting into the nucleus. Sci Signal. 2008 Dec;1(51):pe54.17. Levine JH, Lin Y, Elowitz MB. Functional roles of pulsing in genetic circuits. Science. 2013Dec;342(6163):1193–1200.18. Connor FR. Introductory topics in electronics and telecommunication: modulation. The Univer-sities Press, Belfast; 1973.19. Eldar A, Elowitz MB. Functional roles for noise in genetic circuits. Nature. 2010 Sep;467(7312):167–173.20. Locke JC, Young JW, Fontes M, Hern´andez Jim´enez MJ, Elowitz MB. Stochastic pulse regulationin bacterial stress response. Science. 2011 Oct;334(6054):366–369.21. Yissachar N, Sharar Fischler T, Cohen AA, Reich-Zeliger S, Russ D, Shifrut E, et al. Dynamicresponse diversity of NFAT isoforms in individual living cells. Mol Cell. 2013 Jan;49(2):322–330.22. Lisman JE. Bursts as a unit of neural information: making unreliable synapses reliable. TrendsNeurosci. 1997 Jan;20(1):38–43.23. Mora T, Wingreen NS. Limits of sensing temporal concentration changes by single cells. Phys RevLett. 2010 Jun;104(24):248101–248101.24. Tostevin F, de Ronde W, ten Wolde PR. Reliability of frequency and amplitude decoding in generegulation. Phys Rev Lett. 2012 Mar;108(10):108104–108104.25. Yi TM, Huang Y, Simon MI, Doyle J. Robust perfect adaptation in bacterial chemotaxis throughintegral feedback control. Proc Natl Acad Sci U S A. 2000 Apr;97(9):4649–4653.26. Barkai N, Leibler S. Robustness in simple biochemical networks. Nature. 1997 Jun;387(6636):913–917.27. De Palo G, Facchetti G, Mazzolini M, Menini A, Torre V, Altafini C. Common dynamical featuresof sensory adaptation in photoreceptors and olfactory sensory neurons. Sci Rep. 2013 Feb;3:1251.28. Li Y, Goldbeter A. Pulsatile signaling in intercellular communication. Periodic stimuli are moreefficient than random or chaotic signals in a model based on receptor desensitization. Biophys J.1992 Jan;61(1):161–171.29. Rapp PE, Mees AI, Sparrow CT. Frequency encoded biochemical regulation is more accurate thanamplitude dependent control. J Theor Biol. 1981 Jun;90(4):531–544.830. Tyson JJ, Novak B. Cell cycle: who turns the crank? Curr Biol. 2011 Mar;21(5):185–187.31. McWatters H, Dunlap JC, Millar AJ. Circadian biology: clocks for the real world. Curr Biol. 1999Sep;9(17):633–635.32. Dolmetsch RE, Xu K, Lewis RS. Calcium oscillations increase the efficiency and specificity of geneexpression. Nature. 1998 Apr;392(6679):933–936.33. Batchelor E, Loewer A, Mock C, Lahav G. Stimulus-dependent dynamics of p53 in single cells.Mol Syst Biol. 2011 May;7:488–488.34. Ashall L, Horton CA, Nelson DE, Paszek P, Harper CV, Sillitoe K, et al. Pulsatile stimulation deter-mines timing and specificity of NF- κ B-dependent transcription. Science. 2009 Apr;324(5924):242–246.35. Tay S, Hughey JJ, Lee TK, Lipniacki T, Quake SR, Covert MW. Single-cell NF- κ B dynamics revealdigital activation and analogue information processing. Nature. 2010 Jul;466(7303):267–271.36. Albeck JG, Mills GB, Brugge JS. Frequency-modulated pulses of ERK activity transmit quanti-tative proliferation signals. Mol Cell. 2013 Jan;49(2):249–261.37. Cai H, Katoh-Kurasawa M, Muramoto T, Santhanam B, Long Y, Li L, et al. Nucleocytoplasmicshuttling of a GATA transcription factor functions as a development timer. Science (New York,NY). 2014 Mar;343(6177):1249531.38. Gregor T, Fujimoto K, Masaki N, Sawai S. The onset of collective behavior in social amoebae.Science. 2010 May;328(5981):1021–1025.39. Goldbeter A, G´erard C, Gonze D, Leloup JC, Dupont G. Systems biology of cellular rhythms.FEBS Lett. 2012 Aug;586(18):2955–2965.40. Takeda K, Shao D, Adler M, Charest PG, Loomis WF, Levine H, et al. Incoherent feedforwardcontrol governs adaptation of activated ras in a eukaryotic chemotaxis pathway. Sci Signal. 2012Jan;5(205).41. Shen-Orr SS, Milo R, Mangan S, Alon U. Network motifs in the transcriptional regulation networkof
Escherichia coli . Nat Genet. 2002 May;31(1):64–68.42. Goentoro L, Shoval O, Kirschner MW, Alon U. The incoherent feedforward loop can providefold-change detection in gene regulation. Mol Cell. 2009 Dec;36(5):894–899.43. Hille B. Ion channels of excitable membranes. 3rd ed. Sinauer Associates Inc 2001-07; 2001.44. Endres RG, Wingreen NS. Maximum likelihood and the single receptor. Phys Rev Lett. 2009Oct;103(15):158101–158101.45. Endres RG. Polar chemoreceptor clustering by coupled trimers of dimers. Biophys J. 2009Jan;96(2):453–463.46. Wiggins P, Phillips R. Membrane-protein interactions in mechanosensitive channels. Biophys J.2005 Feb;88(2):880–902.47. Alon U, Surette MG, Barkai N, Leibler S. Robustness in bacterial chemotaxis. Nature. 1999Jan;397(6715):168–171.48. Berg H, Purcell E. Physics of chemoreception. Biophysical Journal. 1977 Nov;20(2):193–219.949. Yoshimoto H, Saltsman K, Gasch AP, Li HX, Ogawa N, Botstein D, et al. Genome-wide analysis ofgene expression regulated by the calcineurin/Crz1p signaling pathway in
Saccharomyces cerevisiae .J Biol Chem. 2002 Aug;277(34):31079–31088.50. Hoffmann A, Baltimore D. Circuitry of nuclear factor κ B signaling. Immunol Rev. 2006Apr;210:171–186.51. Schmitt AP, McEntee K. Msn2p, a zinc finger DNA-binding protein, is the transcriptional acti-vator of the multistress response in Saccharomyces cerevisiae. Proc Natl Acad Sci U S A. 1996Jun;93(12):5777–5782.52. Kang J, Xu B, Yao Y, Lin W, Hennessy C, Fraser P, et al. A dynamical model reveals geneco-localizations in nucleus. PLoS Comput Biol. 2011 Jul;7(7).53. Ferrell JE. Bistability, bifurcations, and Waddington’s epigenetic landscape. Curr Biol. 2012Jun;22(11):458–466.54. Purvis JE, Karhohs KW, Mock C, Batchelor E, Loewer A, G L. p53 dynamics control cell fate.Science. 2012 Jun;336:1440–1444.55. Muramoto T, Cannon D, Gierlinski M, Corrigan A, Barton GJ, Chubb JR. Live imaging of nascentRNA dynamics reveals distinct types of transcriptional pulse regulation. Proc Natl Acad Sci U SA. 2012 May;109(19):7350–7355.56. Thattai M, van Oudenaarden A. Intrinsic noise in gene regulatory networks. Proc Natl Acad SciU S A. 2001 Jul;98(15):8614–8619.57. Lang AH, Fisher CK, Mora T, Mehta P. Thermodynamics of statistical inference by cells. PhysRev Lett. 2014 Oct;113(14):148103–148103.58. Csan´ady L, Vergani P, Gadsby DC. Strict coupling between CFTR’s catalytic cycle and gating ofits Cl- ion pore revealed by distributions of open channel burst durations. Proc Natl Acad Sci US A. 2010 Jan;107(3):1241–1246.59. Schneggenburger R, Ascher P. Coupling of permeation and gating in an NMDA-channel poremutant. Neuron. 1997 Jan;18(1):167–177.60. Arron JR, Winslow MM, Polleri A, Chang CP, Wu H, Gao X, et al. NFAT dysregulation byincreased dosage of DSCR1 and DYRK1A on chromosome 21. Nature. 2006 Jun;441(7093):595–600.61. Anderson PW. More is different. Science. 1972 Aug;177(4044):393–6.62. Mehta P, Schwab DJ. Energetic costs of cellular computation. Proc Natl Acad Sci U S A. 2012Oct;109(44):17978–17982.63. Gillespie DT. A general method for numerically simulating the stochastic time evolution of coupledchemical reactions. J Comput Phys. 1976;22(4):403 – 434.0
Supporting Information Legends
S1 Text.
Details of analytical calculations.
S1 Fig. First three moments of the protein distribution in concentration sensing from themaster equation.
Averages (A,B), variance (C,D), and skewness (E,F) as a function of the frequencyof binding events, f = k + c / (1 + k + c /k − ). (Insets) Magnification of small-noise approximation region(fast switching). Analytical results for CM (blue) and numerical results for BM (red) and intermediatemodulation IM (green) as function of the frequency of binding events (logarithmic scale). Note that thisfigure is similar to Fig. 5 in main text with the addition of IM. Two regimes are shown: k − = 10 k + c ( α = 100 s − , γ = 1 s − , ζ from 1000 to 1) (left column) and k − = 0 . k + c ( α = 10 s − , γ = 1 s − , ζ from 1000 to 1) (right column). Averages from CM, BM and IM are constrained to be equal, i.e. ζ (BM) = αk − − (CM) = α ′ τ b (IM). Variances of CM, BM and IM exhibit two different regimes for fastswitching: for k + c < k − BM is the most accurate mechanism and CM the worst (inset in C) while for k + c > k − CM is generally the most accurate (except for ζ = 1) and IM the worst (inset in D). Thirdmoments show that, for large noise, the probability distributions become asymmetric. S2 Fig. Examples of time traces of receptor activity and protein copy numbers fordifferent regimes. (Top) Regime k + c < k − with k + c = 0 . k − ( α = 100 s − , γ = 1 s − ). (Bottom)Regime k + c > k − with k + c = 10 k − ( α = 10 s − , γ = 1 s − ). (Left) Slow switching with ζ = 400.(Right) Fast switching with ζ = 7. Receptor activity r and protein copy numbers n ( t ) for CM, BM andIM are shown in black, blue, red and green, respectively. S3 Fig. Investigating accuracy based on accumulative signaling (without proteinproduction and degradation). (A) Regime k + c < k − with k + c = 0 . k − ( α = 100 s − , γ = 1 s − and ζ = 7). (Left) ODE model. (Right) Stochastic protein production during τ b in CM and IM. (Top)Examples of time traces. (Bottom) Histograms of number of proteins produced after 100s withstandard deviation in legend based on 1000 simulations. (B) Analogous to (A) but for regime k + c > k − with k + c = 10 k − ( α = 100 s − , γ = 1 s − and ζ = 7). CM, BM and IM are shown in blue, red andgreen, respectively.1 S4 Fig. Incoherent feedforward loop: Comparison of analytical results with simulations ofthe stochastic differential equations. (A) Averages of signaling rate u (left), species y from Eq.(S42) (middle) and species x from (S41) (right) as a function of time. Analytic solutions Eqs. (S32),(S43) and (12) are shown for BM in red, while a (time averaged) time-trace from a stochasticsimulation using the Euler method is shown in orange (CM is almost identical and hence is not shown).(B) Corresponding variances as a function of time for k + c > k − ( k − = 6 . × s − , k + c = 10 s − ).Analytic results are shown in blue for CM and in red for BM; average over time (1 s ) from numericalsimulations are shown in light blue for CM and in orange for BM. (C) Corresponding variances as afunction of time for k + c < k − ( k − = 6 . × s − , k + c = 10 s − ). Colors same as in (B). Remainingparameters: k + c = 10 s − , k x = 10 s − and k y = 50 s − . S5 Fig. Integral feedback loop: Comparison of analytical results with simulations of thestochastic differential equations. (A) Averages of signaling rate u (left), species y from Eq. (S60)(middle) and species x from (S59) (right) as a function of time. Analytic solutions Eqs. (S32), (S66)and (S65) are shown for BM in red, while a (time averaged) time-trace from a stochastic simulationusing the Euler method is shown in orange (CM is almost identical and hence is not shown). (B)Corresponding variances as a function of time for k + c > k − ( k − = 6 . × s − , k + c = 10 s − ).Analytic results are shown in blue for CM and in red for BM; numerical simulations are shown in lightblue for CM and in orange for BM. (C) Corresponding variances as a function of time for k + c < k − ( k − = 6 . × s − , k + c = 10 s − ). Colors same as in (B). Remaining parameters: k + c = 10 s − , k x = 10 s − and k y = 50 s − . S6 Fig. From CM (BM) to AM (FM) for multiple receptors/ion channels. (A-D) Schematicof receptor activity in time. (A) AM emerges from N unsynchronized receptors or ion channels in CMmode. (B) N synchronized CM receptors lead to a hybrid mechanism with information encoded in thefrequency of broad bursts of variable duration. (C) N unsynchronized BM receptors provide a denseseries of bursts. For large N , bursts may start overlapping, leading to variable amplitudes. (D) FMemerges from N synchronized receptors in BM mode. (E) Relative variance for a system of 8 receptorswith ρN synchronized and (1 − ρ ) N unsynchronized receptors, plotted for fast dynamics in the k + c < k − regime (CM in blue and BM in red). Letters refer to panel labels (A-D). Dotted red lineindicates uncertainty from FM for comparison. (Inset) Same for a system of two receptors only. r X i v : . [ q - b i o . S C ] M a y Supporting Informations:
Accurate encoding and decoding by single cells:amplitude versus frequency modulation
Gabriele Micali,
1, 2, 3
Gerardo Aquino,
1, 2
David M. Richards,
1, 2 and Robert G. Endres
1, 2, ∗ Department of Life Sciences, Imperial College, London, UK Centre for Integrative System Biology and Bioinformatics, Imperial College, London, UK Dipartimento di Fisica, Universit`a degli Studi di Milano, Milano, Italy
TEXT S1CONCENTRATION SENSING BY CM RECEPTOR
In this section we calculate the analytic solution for the master equation (Eq. ?? ,b in the main text) for continuousmodulation (CM). For clarity, we repeat here the master equations for the on and off states: dp on ( n, t ) dt = γ ( n + 1) p on ( n + 1 , t ) + αp on ( n − , t ) + k + c p off ( n, t ) − ( γn + α + k − ) p on ( n, t ) ,dp off ( n, t ) dt = γ ( n + 1) p off ( n + 1 , t ) + k − p on ( n, t ) − ( γn + k + c ) p off ( n, t ) , where p on / off ( n, t ) is the probability that the receptor/ion channel is in the on / off state with n output proteins attime t , α and γ are the production and degradation rates respectively, and k + c and k − are the binding and unbindingrates respectively. Note that the concentration of the input species ( c ) is now constant. Eqs. ( ?? ,b) can be rewrittenas Eq. ( ?? ) in the main text dp s ( n, t ) dt = γ ( n + 1) p s ( n + 1 , t ) + α s p s ( n − , t ) + k ¯ s p ¯ s ( n, t ) − ( γn + α s + k s ) p s ( n, t ) , where ¯ s is the on / off state when s is the off / on state, and α on = α , α off = 0, k on = k − and k off = k + c . To find thesolution for the first two moments of the distribution p ( n, t ), we now follow Mehta and Schwab [1]. At steady stateEq. ( ?? ) becomes K ¯ s p ¯ s ( n ) = − ( n + 1) p s ( n + 1) − A s p s ( n −
1) + ( n + A s + K s ) p s ( n ) , (S2)where K s = k s /γ and A s = α s /γ . Using the generating function in Eq. ( ?? ) G s ( z ) = ∞ X n =0 p s ( n ) z n , Eq. (S2) becomes [( z − ∂ z − A s ( z −
1) + K s ] G s ( z ) = K ¯ s , (S3)which implies ( ∂ z − A on ) G on = K off G off − K on G on z − , (S4)( ∂ z − A off ) G off = K on G on − K off G off z − , (S5)which, when combined, gives ( ∂ z − A on ) G on = − ( ∂ z − A off ) G off . (S6) ∗ [email protected] To proceed further, it is useful to define the quantity H s ( z ) related to the generating function G s ( z ) by G s ( z ) = e A s z H s ( z ) . (S7)Using (cid:16) ∂ z − A s (cid:17) G s ( z ) = e A s z ∂ z H s ( z ) , (S8)Eq. (S6) becomes e A on z ∂ z H on ( z ) = − e A off z ∂ z H off ( z ) , (S9)which links the expressions for H on and H off . At this point the initial equation for the steady state (Eq. (S3)) becomes( z − e A s z ∂ z H s ( z ) + K s e A s z ∂ z H s ( z ) = K ¯ s e A ¯ s z ∂ z H ¯ s ( z ) . (S10)Multiplying by e A ¯ s , taking the derivative with respect to z , substituting Eq. (S9), and defining ∆ A s = A ¯ s − A s , gives ∂ z H s ( z ) − ∆ A s ( z − ∂ z H s ( z ) + ( z − ∂ z H s ( z ) − K s ∆ A s H s ( z ) + K s ∂ z H s ( z ) = − K ¯ s H s ( z )and hence ( z − ∂ z H s ( z ) + (cid:16) − ∆ A s ( z −
1) + K s + K ¯ s (cid:17) ∂ z H s ( z ) − K s ∆ A s H s ( z ) = 0 . Finally, changing variables to u = ∆ A s ( z −
1) provides u∂ z H s ( u ) + (1 + K s + K ¯ s − u ) ∂ z H s ( u ) − K s H s ( u ) = 0 . (S11)This is the confluent hypergeometric equation, for which the solution in terms of confluent hypergeometric functionsof the first kind is given by H s ( u ) = c s F ( K s , K s + K ¯ s ; u ) , (S12)with c s a constant of integration. Thus, through Eq. (S7), G s ( z ) = c s e A s z F ( K s , K s + K ¯ s ; ∆ A s ( z − . (S13)To determine the constants, notice that F ( a, b,
0) = 1 leading to G s (1) = c s e A s = h p s i = K ¯ s K ¯ s + K s , where h p s i is the average probability of being in state s . Rearranging terms, we obtain c s = K ¯ s K ¯ s + K s e − A s . (S14)Finally, the probability distribution at steady state is given by [1] G s ( z ) = K ¯ s e A s ( z − K ¯ s + K s F ( K s , K s + K ¯ s ; ∆ K s ( z − . (S15)Having an analytic expression for the steady-state probability distribution (Eq. S15), we can now calculate the first,second and third moments, which are related to the mean, variance and skewness, respectively. The mean productionof the output protein is given by the mean production in the on state multiplied by the probability to be in the on state, averaged over the whole time period. For such a two-state system h p on i = K off K off + K on and h p off i = 1 − h p on i .Therefore, the mean number of proteins is given by h n i = ( A on − A off ) h p on i + A off = αγ k + c/k − k + c/k − . (S16)To calculate the variance, we use the following property of the generating function:( δn ) = X s ( ∂ z z∂ z G s ( z )) (cid:12)(cid:12) z =1 − h n i . (S17) Proof. ( δn ) = X s "X n n p s ( n ) − h n i = X s "X n n p s ( n ) z =1 − h n i = X s "X n z∂ z ( np s ( n ) z n ) n p s ( n ) z =1 − h n i = X s [ z∂ z z∂ z G s ( z )] z =1 − h n i = X s [ ∂ z z∂ z G s ( z )] z =1 − h n i . Using common properties of hypergeometric functions, the analytical solution for the variance is [1] D ( δn ) E = h n i + h p on i h p off i (∆ A s ) K s + K ¯ s = h n i + α γ ( γ + k − + k + c ) k + c/k − (1 + k + c/k − ) . (S18)For details, see the full calculation in the SI of [1]. Third moment
In order to understand more about the symmetry of the probability distribution, we calculate the third moment atsteady state. As in Eq. (S17) the third moment can be found via generating functions as (cid:10) n (cid:11) = X s [ ∂ z z∂ z z∂ z G s ( z )] (cid:12)(cid:12)(cid:12) z =1 = X s (cid:2) z∂ z G s ( z ) + z ∂ z G s ( z ) + ∂ z G s ( z ) (cid:3) (cid:12)(cid:12)(cid:12) z =1 , (S19)where [1] X s ∂ z G s ( z ) (cid:12)(cid:12)(cid:12) z =1 = h n i , (S20) X s (cid:0) z∂ z G s ( z ) (cid:1) (cid:12)(cid:12)(cid:12) z =1 = h p on i ( A on ) + h p off i ( A off ) − h p on i h p off i ∆ A s ( K on + K off )1 + K s + K ¯ s = α ( k + c/k − ) γ (1 + k + c/k − ) + α ( k + c ) γ (1 + k + c/k − )( γ + k − + k + c ) . (S21)Thus, only P s z ∂ z G s ( z ) (cid:12)(cid:12) z =1 needs to be calculated. The result is X s (cid:0) z ∂ z G s ( z ) (cid:1) (cid:12)(cid:12)(cid:12) z =1 = ( A on ) h p on i + ( A off ) h p off i − h p on i h p off i (∆ A on ) ( K on + K off )1 + K on + K off × (cid:20) (3 + K on + 2 K off ) A on K on + K off + (3 + 2 K on + K off ) A off K on + K off (cid:21) = α ( k + c/k − ) γ (1 + k + c/k − ) − α ( k + c/k − ) (3 γ + k − + 2 k + c ) γ (1 + k + c/k − ) ( γ + k − + k + c ) (2 γ + k − + k + c ) . (S22)By combining Eqs. (S57)-(S59) as indicated in Eq. (S56), we obtain the analytic expression for the skewness of oursystem. Proof.
From the definition of confluent hypergeometric functions of the first kind F ( a, b ; z ) = ∞ X n =0 a n b n n ! z n , with a n = a ( a + 1)( a + 2) ... ( a + n − ,b n = b ( b + 1)( b + 2) ... ( b + n − , and the fact that z∂ z F ( a, b ; z ) = z ab F ( a + 1 , b + 1; z ) , (S23)we obtain ∂ z G s ( z ) = K ¯ s K ¯ s + K s h A s e A s ( z −
1) 1 F ( · , · , ∆ A s ( z − A s K s K s + K ¯ s e A s ( z −
1) 1 F (+ , + , ∆ A s ( z − i , (S24)where ∆ A s = A ¯ s − A s and F (+ , + , z ) = F ( a + 1 , b + 1 , z ). We now need to calculate P s z ∂ z G s ( z ) (cid:12)(cid:12) z =1 . UsingEq. (S21) we find that (cid:0) z ∂ z G s ( z ) (cid:1) (cid:12)(cid:12)(cid:12) z =1 = z K ¯ s K ¯ s + K s ∂ z (cid:20)(cid:16) ( A s ) F ( · , · , ∆ A s ( z − A s K s A s K s + K ¯ s F (+ , + , ∆ A s ( z − A s ) K s ( K s + 1)(1 + K s + K ¯ s )(2 + K s + K ¯ s ) · F (++ , ++ , ∆ A s ( z − (cid:17) e A s ( z − (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) z =1 . (S25)Differentiating and using property (S23), we obtain (cid:0) z ∂ z G s ( z ) (cid:1) (cid:12)(cid:12)(cid:12) z =1 = z K ¯ s K ¯ s + K s e A s ( z − h ( A s ) F ( · , · , ∆ A s ( z − A s K s ( A s ) K s + K ¯ s F (+ , + , ∆ A s ( z − A s ) K s ( K s + 1) K s (1 + K s + K ¯ s )(2 + K s + K ¯ s ) F (++ , ++ , ∆ A s ( z − A s ) K s ( K s + 1) ( K s + 2)(1 + K s + K ¯ s )(2 + K s + K ¯ s )(3 + K s + K ¯ s ) F (+ + + , + + + , ∆ A s ( z − i(cid:12)(cid:12)(cid:12)(cid:12) z =1 . (S26)Evaluating at z = 1, we obtain (cid:0) z ∂ z G s ( z ) (cid:1) (cid:12)(cid:12)(cid:12) z =1 = K ¯ s K ¯ s + K s " ( A s ) + ∆ A s K s K s + K ¯ s (cid:18) A s ) + ∆ A s ( K s + 1)2 + K s + K ¯ s (cid:18) A s + ∆ A s ( K s + 2)3 + K s + K ¯ s (cid:19) (cid:19) . (S27)Finally, summing on the possible states of s we arrive at Eq. (S22). SMALL-NOISE APPROXIMATION TO RAMP SENSINGInput noise
In the Model section of the main text, we built a model for a single receptor/ion channel that encodes informationfrom an cell-external environment in some cell-internal degrees of freedom. Similarly to [2], we assume that thereceptor/ion channel activity ( r ( t )) is a two state system: on with r = 1 when the receptor is bound or the channelopen, and off with r = 0 when the receptor is not bound or the channel is closed. The external concentration( c ( t )) is assumed to affect the unbound/closed time interval h τ u i = [ k + c ( t )] − but not the bound/open time internal h τ b i = k − − , where k + and k − are both constants. Both interval durations are assumed to be independent, exponentiallydistributed random variables. The independence of binding and unbinding (or equivalently of opening and closing)means that the probability of a molecule binding the receptor a second time is negligible. We therefore assume thesystem to be in the fast diffusion regime.The signaling rate, called u , implements two different mechanisms of encoding, either continuous (CM) or bursty(BM) modulation. CM and BM ultimately correspond to amplitude (AM) and frequency (FM) modulation, respec-tively, when generalized to multiple receptors/ion channels as explained in the Results section of the main text. InCM the proteins are produced with a constant rate α during the binding time. On the other hand, for BM a burst of ζ proteins is realized at the time of binding, so u ( t ) = ( α r ( t ) for CM ζ P δ ( t − t + i ) for BM, (S28)where αk − − = ζ and t + i the binding times. By taking the average of the rate u ( t ) over a time ˜ t much longer thanboth the average bound time, h τ b i = k − − , and the average unbound time, h τ u i = [ k + c ( t )] − , but shorter than thetime during which the external concentration changes, we obtain h u ( t ) i = ζ k + c ( t )1 + k + c ( t ) /k − , (S29) h δu ( t ) i = 0 , h δu ( t ) δu ( t ′ ) i = gζ k + c ( t )(1 + k + c ( t ) /k − ) δ ( t − t ′ ) , (S30) g = ( k + c ( t ) /k − ] for BM , (S31)which become Eqs. ( ?? ), ( ?? ) and ( ?? ) in the main text by setting ζ = 1. (See Supplementary Information in [2] forfurther details). Importantly, g BM < g CM since k + c ( t ) < k − , i.e. h τ b i < h τ u i .By considering an external concentration given by Eq. ( ?? ) in the main text, c ( t ) = ( c + c t, t ≥ c , t < , with c t ≪ c , Eq. (S29) becomes u ( t ) = ζ ( u + u t + δu ( t ) , t ≥ u + δu ( t ) , t < , (S32)where we assume δu ( t ) ≪ u + u t and u = k + c (1 + k + c /k − ) , (S33) u = k + c (1 + k + c /k − ) , (S34) h δu ( t ) δu ( t ′ ) i t,t ′ ≥ = ( δ ( t − t ′ )(1+ k + c /k − ) [ g CM u + g CM (1 − k + c /k − ) u t ] for CM δ ( t − t ′ )(1+ k + c /k − ) (cid:2) g ∗ BM u + (cid:0) − k + c /k − + 3 k c /k − (cid:1) u t (cid:3) for BM , (S35) h δu ( t ) δu ( t ′ ) i t,t ′ < = gu (1 + k + c /k − ) δ ( t − t ′ ) . (S36)Here g CM = 2 and g ∗ BM = 1 + ( k + c /k − ) (cf. Eq. ( ?? ) in the main text). Again for ζ = 1, this becomes Eq. ( ?? )in the main text. From Eqs. (S32)-(S36), the constant ( t <
0) and ramp ( t ≥
0) regimes for the external species c are encoded in the rate u in the corresponding regimes since the condition c t ≪ c ensures u t ≪ u . However, tosatisfy condition δu ≪ h u i for both CM and BM, a new condition is needed: g CM / BM k + c (1 + k + c /k − ) ≪ τ c , (S37)which implies k − , k + c ≫ τ − c . (S38)Here, we have introduced the correlation time of white noise, τ c , corresponding to the δ -function used in Eqs.(S35)and (S36). Note that condition in (S38) restricts our study to the fast switching regime. Finally, the signaling rate u in the constant regime has one small term δu/u of order δ which is defined by Eqs. (S32) and (S36). Instead,small-ramp regime u contains two small terms: the small-ramp term u t/u of order ǫ and the small noise term δu/u which now has a correction to order δ coming from the small ramp (order ǫ ). With these definitions, from Eqs.(S32)-(S35) the rate u in the small ramp regime has two small corrections to the constant rate u u ( t ) = ζu u tu |{z} o ( ǫ ) + δuu |{z} o ( δ (1+ √ ǫ ) ) ! , (S39)In order to linearize around linear solutions, we further assume that the small-noise amplitude is smaller than thesmall ramp. As a result, o ( δ ) ∼ o ( ǫ x ) with x >
1, which means that g (1 + k + c /k − ) k + c < τ c ( k + c t ) . (S40)Note that for simplicity, both in the following sections and in the main text, we set ζ = 1. Output noise in incoherent feedforward loop
Average solutions for ramp sensing
Eqs. ( ?? ) and ( ?? ) in the main text for the incoherent feedforward loop for f ( u ) = e bu and g ( y ) = e − bk y y become dxdt = k x h e b ( u − k y y ) − x i , (S41) dydt = u − k y y. (S42)Here, b is a constant introduced to maintain the exponent unitless, and k x and k y are rate constants for x and y .This system of equations performs exact adaptation. The steady-state solution in the constant regime ( t < h x ( t ) i = 1 , h y ( t ) i = u k y , (S43)which sets the initial conditions h x (0) i = 1 and h y (0) i = u /k y for Eqs. (S41) and (S42) in the ramp regime ( t ≥ ?? ). With these initial conditions, the solutions for t ≥ h x ( t ) i = e − k x t k x e bu ky Z t dt ′ e k x t ′ − bu ky e − kyt ′ , (S44) h y ( t ) i = u k y − u k y + u k y t + u k y e − k y t , (S45)where the integral in the expression for h x ( t ) i cannot be solved analytically. However, by assuming that the integralstarts from time ǫ ≫ k − y , Eq. (S44) becomes h x ( t ) i ≃ e bu ky + e − k x ( t − ǫ ) (cid:20) h x ( ǫ ) i − e bu ky (cid:21) . Finally, by considering t such that t − ǫ ≫ k − x and without exceeding the small-ramp regime ( e.g. k x,y ≫ ǫ small), the solution becomes Eqs. ( ?? -b) in the main text, h x ( t ) i = e bu ky , h y ( t ) i = u k y − u k y + u k y t. These solutions match numerical results shown in S4A Fig. Note that the time interval over which these solutionsare valid extends from a time larger than the transient time to around a time that does not exceed the small-rampregime. These criteria also set the regime of validity for our next results.
Output variances in ramp sensing
Above we gave the average solutions for the incoherent feedforward loop, both in the constant regime ( t <
0, Eq.(S43)) and in the ramp regime ( t ≥
0, Eqs. (S44) and (S45)). Now we want to linearize the equations around thesesolutions in order to obtain information about the noise. We assume that the input noise ( δu ) is smaller than theramp, Eqs. (S38) and (S40), which translates into small output noise ( δx and δy ). In addition to these assumptions,we also assume b ∼ u − in order to ensure b ( δu − k y δy ) ≪
1. Hence, in the ramp regime, the differential equationsfor δx and δy become d ( δx ) dt = k x (cid:20) be u bky ( δu − k y δy ) − δx (cid:21) , (S47) d ( δy ) dt = δu − k y δy. (S48)By defining X ( t ) = (cid:20) δxδy (cid:21) , Eqs. (S47) and (S48) can be rewritten in a compact way for both the constant andsmall-ramp regimes as dX ( t ) dt + M ( t ) X ( t ) = (cid:20) A ( t ) δu ( t ) δu ( t ) (cid:21) , (S49)where M ( t ) = " k x k x k y b k y , t < " k x k x k y b e bu ky k y , t ≥ , (S50)and A ( t ) = ( k x b, t < k x b e bu ky , t ≥ . (S51)Note that M ( t ) (for t >
0) is independent of time, which allows Eq. (S49) to be solved analytically. This is due toour choice of f ( u ) and g ( u ) in Eqs. ( ?? ) and ( ?? ) in the main text.For the constant input regime, t <
0, the solutions for CM and BM are D ( δx ( t )) E = g CM/BM b k x u k x + k y )(1 + k + c /k − ) , (S52) D ( δy ( t )) E = g CM/BM u k y (1 + k + c /k − ) , (S53)where g CM = 2 and g BM = g ∗ BM = 1 + ( k + c /k − ) . Hence, in the constant regime the output noise for BM is lowerthan the output noise for CM since k + c < k − .For the small-ramp regime, t ≥
0, Eq. (S49) is analytically solvable for t ≫ k − y by evaluating the integral from time ǫ ≫ k − y to some time t that does not exceed the small-ramp approximation (as discussed for the average solutions).With these assumptions and by using an appropriate integrating factor, the solution for X ( t ) is Z tǫ e Mt ′ dX ( t ′ ) dt ′ dt ′ + Z tǫ e Mt ′ M X ( t ′ ) dt ′ = Z tǫ e Mt ′ (cid:20) Aδu ( t ′ ) δu ( t ′ ) (cid:21) dt ′ ,X ( t ) = e − Mt Z tǫ e Mt ′ (cid:20) Aδu ( t ′ ) δu ( t ′ ) (cid:21) dt ′ + e − M ( t − ǫ ) X ( ǫ ) . However, for t − ǫ ≫ k − x,y (within the limit for t and ǫ as discussed for Eqs. ( ?? ,b) the solution is X ( t ) = e − Mt Z tǫ e Mt ′ (cid:20) Aδu ( t ′ ) δu ( t ′ ) (cid:21) dt ′ . (S54)By using matrix diagonalization, expressing the noise in u by a delta function in time (Eq. S35), and integrating Eq.(S54) for X ( t ) , we find analytical solutions for the variances. The results for CM are D ( δx ( t )) E CM = b k x e bu ky k x + k y )(1 + k + c /k − ) (cid:20) g CM u − g CM (1 − k + c /k − ) k x + k y u + g CM (1 − k + c /k − ) u t (cid:21) , (S55) D ( δy ( t )) E CM = 12 k y (1 + k + c /k − ) (cid:20) g CM u − g CM (1 − k + c /k − ) ky u + g CM (1 − k + c /k − ) u t (cid:21) , (S56)where g CM = 2. For BM, the g BM parameter (see Eq. ( ?? )) affects the integration, and the results are D ( δx ( t )) E BM = b k x e bu ky k x + k y )(1 + k + c /k − ) (cid:20) g ∗ BM u − − k + c /k − + 3 k c /k − k x + k y u + (cid:0) − k + c /k − + 3 k c /k − (cid:1) u t (cid:21) , (S57) D ( δy ( t )) E BM = 12 k y (1 + k + c /k − ) (cid:20) g ∗ BM u − − k + c /k − + 3 k c /k − ky u + (cid:0) − k + c /k − + 3 k c /k − (cid:1) u t (cid:21) , (S58)where g ∗ BM = 1 + ( k + c /k − ) . Note that for c = u = 0, the solutions coincide with the solutions for the constantregime. Furthermore, by comparing the time-dependent terms, CM is noisier than BM when k + c < k − . In ourmodel this is due to the input noise (cf. Eq. S35). These two regimes in which BM is less noisy and hence moreaccurate than CM depend on the ratio of binding and unbinding rates as shown in S4B,C Fig. Clearly the analyticsolutions match numerical simulations with noise. All these calculations were done using Wolfram Mathematica 8,while all the simulations of the stochastic differential equations (Eqs. S41 and S42) were done using the Euler methodin MATLAB. Output noise for integral feedback loop
A similar approach can be applied to the integral feedback loop given by Eqs. ( ?? ) and ( ?? ) in the main text,shown here for clarity with f ( y ) = e − by : dxdt = ue − by − k x x, (S59) dydt = k y ( x − . (S60)Assuming that u is given by Eq. (S32), this system does not have analytic solutions in the ramp regime. However, inthe small-ramp regime it is possible to linearize around the solutions of the constant regime. Hence, h x ( t ) i = x + ǫ x ( t )and h y ( t ) i = y + ǫ y ( t ), where x = 1 and y = ln ( u /k x ) b are the solutions for the constant regime with h u i = u .Note that the condition k x < u is required. By linearization, Eqs. ( ?? ) and ( ?? ) become dǫ x ( t ) dt = k x u tu − k x ǫ x ( t ) − k x bǫ y ( t ) − bk x u u ǫ y ( t ) t, (S61) dǫ y ( t ) dt = k y ǫ x ( t ) . (S62)Combining both equations and neglecting the second-order term u t/u ǫ y ( t ), it possible to find a second-order differ-ential equation for ǫ x ( t ), given by d ǫ x ( t ) dt + k x dǫ x ( t ) dt + bk x k y ǫ x ( t ) = k x u u . (S63)The solution is ǫ x ( t ) = u bk y u + C exp " − k x − r k x − bk x k y ! t + C exp " − k x r k x − bk x k y ! t , (S64)with ǫ x ( t ) → u bk y u after a transient time defined by the exponential terms for any k x / − bk x k y . Furthermore, thereare two integration constants C and C . From Eq. (S62) we obtain ǫ y = u bu t − u b k y u . Finally, the solutions oflinearized Eqs. ( ?? ) and ( ?? ) in the small-ramp regime after the transient time are [2] h x ( t ) i = 1 + u bk y u , (S65) h y ( t ) i = y − u b k y u + u bu t. (S66)Within the small-noise approximation (Eq. S38), we want to find expressions for the variances. In the constantregime, u ( t ) = u + δu ( t ) ( t < x ( t ) = 1 + δx ( t ) and y ( t ) = y + δy ( t ). Therefore, the equations forthe noise terms become d ( δx ) dt = k x (cid:20) δuu − bu δy − δx (cid:21) , (S67) d ( δy ) dt = k x δx. (S68)Proceeding similarly to the incoherent feedforward loop, in the constant regime the solution for the variances are [2] D ( δx ) E = g CM/BM k x k + c /k − ) u , (S69) D ( δy ) E = g CM/BM k y b (1 + k + c /k − ) u , (S70)with g CM = 2 and g BM = g ∗ BM = 1 + ( k + c /k − ) . Hence, in the constant regime, the output noise for BM is lowerthan the output noise for CM (since k + c < k − ).To study the system in the small-ramp regime ( t > δu ) is smaller thanthe ramp (Eqs. S38 and S40), which translates into small output noise ( δx and δy ), and linearize around solutions(S65) and (S66). As a result, Eq. ( ?? ) becomes d ( δx ) dt = k x (cid:20)(cid:18) u tu + δuu (cid:19) (cid:18) − b ( ǫ y + δy ) + b ǫ y + δy ) (cid:19) − − ǫ x − δx (cid:21) , which, by using Eq. (S61), becomes d ( δx ) dt = k x b " δuu − δy − u u ǫ y t − u tu δy − ǫ y u δu + bǫ y bǫ y δy − δxb + o ( ǫ , δ , ǫ δ ) , (S71)where we neglect third-order terms in the small ramp ( o ( ǫ )), second-order terms in the small noise ( δ ) and mixed-order terms ( o ( ǫ δ )) due to the assumption that the noise is smaller than the ramp (cf. discussion that leads to Eq.(S40)).Defining X ( t ) = (cid:20) δxδy (cid:21) , Eqs. (S71) and (S60) become dX ( t ) dt + M ( t ) X ( t ) = (cid:20) A ( t ) δu ( t ) + B ( t )0 (cid:21) , (S72)where from Eqs. (S71) and (S60), using definitions of ǫ x and ǫ y , M ( t ) = " k x k x b (cid:16) bu k y u (cid:17) − k y , A ( t ) = k x u (cid:16) u k y bu − u u t (cid:17) and B ( t ) = k x u u (cid:16) tk y − t (cid:17) . Using an integrating factor and integrating between ǫ and t gives X ( t ) = e − M ( t − ǫ ) X ( ǫ ) + e − Mt Z tǫ dz e Mz A ( z ) (cid:20) δu (cid:21) + e − Mt Z tǫ dz e Mz B ( z ) (cid:20) (cid:21) , (S73)0where the term e − M ( t − ǫ ) X ( ǫ ) is negligible for ( t − ǫ ) ≫ k − x,y . To calculate the variances we square Eq. (S73). UsingEqs. (S5)-(S9), the results for D ( δx ( t )) E to first-order in the small-ramp parameters are D ( δx ( t )) E CM = kx k + c /k − ) u (cid:20) g CM − g CM (1 + 2 k + c /k − ) u tu + g CM k x + bk y (1 + 2 k + c /k − ) u bk y k x u (cid:21) , (S74) D ( δx ( t )) E BM = kx k + c /k − ) u (cid:20) g ∗ BM − (cid:0) k + c /k − − k c /k − − (cid:1) u tu + (cid:0) bk y (cid:0) k + c /k − − k c /k − − (cid:1) + 2 k x g ∗ BM (cid:1) u bk x k y u (cid:21) . (S75)Similarly the results for D ( δy ) E are D ( δy ( t )) E CM = k y b (1 + h τ b i k + c ) u " g CM − g CM (1 + 2 k + c /k − ) u tu + g CM ( kx (3 + 2 k + c /k − ) + 2 k y b (1 + 2 k + c /k − )) u bk x k y u , (S76) D ( δy ( t )) E BM = k y b (1 + h τ b i k + c ) u " g ∗ BM − (cid:16) h τ b i k + c − h τ b i k c (cid:17) u tu + h bk y (cid:16) h τ b i k + c − h τ b i k c (cid:17) + k x (cid:16) h τ b i k + c + h τ b i k c (cid:17)i u bk y k x u , (S77)where g ∗ BM = 1 + ( k + c /k − ) . Note that for c = u = 0, the solutions coincide with the solutions for the constantregime. Although it is clear that BM is less noisy than CM for k + c < k − , by comparing the time-dependent termswe find that, in fact, BM is always less noisy than CM. The analytical solutions are plotted in S5 Fig. and matchthe numerical simulations with noise. Again, all these calculations were done using Wolfram Mathematica 8, whileall the simulations of the stochastic differential equations (Eqs. S59 and S60) were done using the Euler method inMATLAB. FURTHER INVESTIGATIONS INTO THE ACCURACY
In this section we provide further explanations for the accuracy of concentration sensing by a single receptorwithout comparing with the maxmimum-likelihood estimation [3]. In Fig. ?? we showed results from the masterequation for the two regimes k + c < k − and k + c > k − for slow and fast switching of the receptor. Despite itsburstiness, the BM receptor turned out more accurate than the CM receptor in the k + c < k − regime for fast switching. Additional results from the master-equation model.
To understand this result better we also implemented anintermediate-modulation (IM) receptor, which has features of both the CM and BM receptors. Like the CM receptor,the IM receptor signals while in the bound ( on ) state, but instead of a constant rate α of production it producesprotein with a rate α ′ so that in each bound interval the same number of molecules are produced irrespective of theinterval length, i.e. α ′ τ b = ζ , with ζ the constant burst size of BM. For this to work, the IM receptor would haveto know at the time of binding when it will unbind again, in order to choose the correct rate of production. Sincethe rate of unbinding is a random variable this is generally not possible. Nevertheless, the IM receptor may help tofurther elucidate our observed trends in accuracy. In practice, we implemented this IM receptor by first simulatinga time trace of bound and unbound time intervals with a Gillespie algorithm, allowing us to determine the rate ofproduction as a function of time. Afterwards, the actual protein production and degradation were simulated.In analogy to Fig. ?? the results for the IM receptor are shown in S1 Fig. (green lines), which also shows theresults for the CM and BM receptors for comparison in blue and red, respectively. As expected, for slow switching theIM receptor has intermediate accuracy between CM and BM. CM is most accurate as continuous production duringthe bound intervals is balanced by degradation so the output protein level does not fluctuate excessively. BM is leastaccurate due to the increased burst size for slow switching. Since signaling by the IM receptor is only burst-like forthe short bound intervals but not for the long bound intervals, it is somewhat more accurate than BM. Due to the1non-constant rate of production, IM also fluctuates more than CM. This intermediate accuracy is clearly demonstratedby the time traces in the left panels of S2 Fig.In the k + c < k − regime for fast switching, the inset of S1C Fig. shows that BM is now most accurate and that IMhas again intermediate accuracy. While BM steadily produces the same amount of protein at the times of binding,IM produces this amount only during short bound intervals as its rate of production is then high, while during longbound intervals its slow production is buffered by degradation, so its protein level fluctuates more strongly. CM iseven worse than IM since, due to its constant rate of production during bound time intervals, it hardly produces anyprotein during short bound intervals, which leads to drastic drops in protein level, while it produces a lot during longbound intervals due to its constant rate of production.In contrast, in the k + c > k − regime for fast switching, CM is generally most accurate due to its approximatelyconstant rate of production throughout time, i.e. the receptor is almost always bound and active. IM is less accuratethan CM because its rate of protein production is variable due to the variable length in bound intervals, despitethe fact that the receptor is mostly bound. Interestingly, IM is even less accurate than BM under these conditions.Inspecting the examples of time trace in the bottom right panel of S2 Fig., the burst sizes of IM can exceed the burstsizes of BM for unusually short bound intervals since production is very high and stochastic, and only on average thesame amount of protein is produced during bound intervals than during a burst in BM. During long bound intervalsthe rate of production is very low. Hence, compared to BM, degradation prevents a net increase in protein levelduring a bound interval, leading to further variability. A special case is when the burst size ζ is 1. As shown in theinset of Fig. ?? D, BM can be more accurate than CM. This is because the burst size of BM is minimal and in themaster equation the production with minimal rate α in CM is highly stochastic.As we now discuss, to provide further intuition for the differences in accuracy between the k + c < k − and k + c > k − regimes, we also simulated the variance of the signaling output (and hence the accuracy-determining factor g ) directly(see Eq. 5). Signaling output from ODE model without protein production and degradation.
Factor g in Eq. 5 (andEq. S31) determines the variance of the signaling rate u ( t ) without invoking any downstream protein productionand degradation. For a given time interval ∆ t , we can hence simulate u ( t ) directly. We assess the accuracy of CM,IM, and BM by plotting the histograms of the integrated signaling rate u I (∆ t ) := R ∆ t u ( t )d t and by determing theirvariances (cf. derivation of g in [2]). As slow protein production and degradation strongly affect the accuracy of thefinal protein output for slow switching, this approach mainly helps understand the interesting fast switching case.We initially assume signaling during bound intervals is deterministic, leading to a linear increase of u with slope α ( α ′ ) during a bound time interval for CM (IM) and a step increase by ζ for BM. At each unbound time interval, IMand BM have the same level of signaling output as IM produces the same number of proteins deterministically duringeach bound interval ( ζ ). In contrast, the signaling output from CM is generally different since the rate of signaling isalways the same for each bound interval but their durations vary. Resulting time traces and variances are shown inS3A and B Figs. left panels, respectively. Specifically, S3A Fig., left panels shows clearly that for k + c < k − BM andIM are most accurate with u I ( t ) increasing almost linearly in time. Since signaling is deterministic, BM and IM areessentially identical, and their variance may only differ due to small differences in signaling during the final boundinterval (S3A Fig., bottom left panel). This last bound time interval may be interrupted in IM, but for long ∆ t thisdifference is negligible. In contrast, S3B Fig., left panels show clearly that for k + c > k − CM is most accurate, as u I ( t ) is now almost linear in time. Signaling output from master-equation model without protein production and degradation.
Allowingsignaling to be stochastic does not change the results for the accuracy significantly. S3A Fig. right panels show thatfor k + c < k − BM is now most accurate and that IM has intermediate accuracy (between BM and CM) due to itsvariability in signaling in line with S1C Fig. Additionally, S3B Fig., right panels show that CM is still most accuratebut also that IM is worse than BM in line with S1D Fig.Taken together, these additional simulation results confirm our findings of the main text that BM is most accuratefor k + c < k − and CM is generally most accurate for k + c > k − . AM IS MORE ACCURATE THAN FM FOR MULTIPLE RECEPTORS/ION CHANNELS
Here, we provide a more detailed discussion of the accuracy of encoding by multiple receptors, i.e. using AM andFM. To determine whether AM or FM is more accurate in encoding and decoding, we generalize to multiple receptors(or ion channels) (S6 Fig.). We assume that AM is obtained by unsynchronized CM receptors (S6A Fig.), while FMis obtained by synchronized receptors that individually operate with BM (S6D Fig.). Other types of synchronizationare also possible with synchronized CM receptors shown in S6B Fig. and unsynchronized BM receptors shown in2S6C Fig. However, these receptors exhibit imperfect FM: resulting bursts have either variable duration (S6B Fig.) orvariable amplitude (S6C Fig.) in contrast to the data (Fig. ?? ) [4, 5].To estimate the accuracy, we first consider perfect synchronization and unsynchronization in either modulationscheme. For N unsynchronized ( us ) receptors, we can express the resulting average and variance of the encodedinput by the single-receptor quantities, i.e. h u ( t ) i usN = N h u ( t ) i and h δu ( t ) δu ( t ′ ) i usN = N h δu ( t ) δu ( t ′ ) i . As a result,the relative variance (variance divided by the average-squared) scales with N − . In contrast, for N synchronized ( s )receptors, the average and variance of the encoded input can be written as h u ( t ) i sN = N h u ( t ) i and h δu ( t ) δu ( t ′ ) i sN = N h δu ( t ) δu ( t ′ ) i , respectively. Hence, the relative variance is now independent of N , so unsynchronized receptorshave an N times smaller noise than synchronized receptors. Since N unsynchronized CM receptors lead to AM, weobtain for its relative variance h δu ( t ) δu ( t ′ ) i AMN (cid:16) h u ( t ) i AMN (cid:17) = h δu ( t ) δu ( t ′ ) i CM N (cid:16) h u ( t ) i CM (cid:17) . (S78)Conversely since N synchronized BM receptors lead to FM, the relative variance of FM is h δu ( t ) δu ( t ′ ) i F MN (cid:16) h u ( t ) i F MN (cid:17) = h δu ( t ) δu ( t ′ ) i BM (cid:16) h u ( t ) i BM (cid:17) . (S79)For slow dynamics, or fast dynamics with k + c > k − , CM is more accurate than BM. Hence, for N receptors, AM iseven more accurate than FM. In contrast, for fast dynamics with k + c < k − , BM is up to twice as accurate as CM(Eq. ( ?? )), and AM is N times more accurate than CM. Consequently, AM becomes more accurate for encoding thanFM for more than two receptors (S6E Fig.). An exception are two receptors, for which AM and FM can be equallyaccurate (S6E Fig., inset). Since we generally show that larger signaling noise leads to larger output noise, the samerule emerges for decoding.To extend our results to intermediate levels of synchronization for N > ρ ofsynchronized receptors while the remaining fraction (1 − ρ ) are unsynchronized, with signaling either by CM or BM(S6E Fig.). When comparing CM and BM receptors for the same levels of synchronization ρ , BM receptors canremain more accurate than CM receptors (S6E Fig.). However, intermediate levels of synchronization do not strictlyrepresent AM and FM. As shown in S6B,C Figs. synchronized CM receptors lead to pulses of variable duration, whileunsynchronized BM receptors lead to highly frequent pulses with potentially variable amplitude.Taken together, since single cells have thousands of receptors and ion channels, AM is the most accurate modulationscheme.3 [1] Mehta P, Schwab DJ (2012) Energetic costs of cellular computation. Proc Natl Acad Sci U S A 109: 17978-17982.[2] Mora T, Wingreen NS (2010) Limits of sensing temporal concentration changes by single cells. Phys Rev Lett 104: 248101-248101.[3] Endres RG, Wingreen NS (2009) Maximum likelihood and the single receptor. Phys Rev Lett 103: 158101-158101.[4] Cai L, Dalal CK, Elowitz MB (2008) Frequency-modulated nuclear localization bursts coordinate gene regulation. Nature455: 485-490.[5] Hao N, O’Shea EK (2012) Signal-dependent dynamics of transcription factor translocation controls gene expression. NatStruct Mol Biol 19: 31-39. h n i + ( δ n ) , Fr e q u e n c y f + ( δ n ) , Fr e q u e n c y f E FCA BD
Slow FastFastSlow k + c < k − k + c > k − S1 Fig.
First three moments of the protein distribution in concentration sensing from the master equation.
Averages (A,B), variance (C,D), and skewness (E,F) as a function of the frequency of binding events, f = k + c / (1 + k + c /k − ).(Insets) Magnification of small-noise approximation region (fast switching). Analytical results for CM (blue) and numericalresults for BM (red) and intermediate modulation IM (green) as function of the frequency of binding events (logarithmic scale).Note that this figure is similar to Fig. ?? in main text with the addition of IM. Two regimes are shown: k − = 10 k + c ( α = 100 s − , γ = 1 s − , ζ from 1000 to 1) (left column) and k − = 0 . k + c ( α = 10 s − , γ = 1 s − , ζ from 1000 to 1) (rightcolumn). Averages from CM, BM and IM are constrained to be equal, i.e. ζ (BM) = αk − − (CM) = α ′ τ b (IM). Variances ofCM, BM and IM exhibit two different regimes for fast switching: for k + c < k − BM is the most accurate mechanism and CMthe worst (inset in C) while for k + c > k − CM is generally the most accurate (except for ζ = 1) and IM the worst (inset inD). Third moments show that, for large noise, the probability distributions become asymmetric. n ( t ) Time [s] n ( t ) Time [s] r CMIMBM
Slow Fast
01 01 r ( t ) r ( t ) k + c < k − k + c > k − S2 Fig.
Examples of time traces of receptor activity and protein copy numbers for different regimes. (Top) Regime k + c < k − with k + c = 0 . k − ( α = 100 s − , γ = 1 s − ). (Bottom) Regime k + c > k − with k + c = 10 k − ( α = 10 s − , γ = 1 s − ).(Left) Slow switching with ζ = 400. (Right) Fast switching with ζ = 7. Receptor activity r and protein copy numbers n ( t ) forCM, BM and IM are shown in black, blue, red and green, respectively. Time [s] u ( t ) CMIMBM500 1000 15000100200300 u ( t =100s) σ ≈ σ ≈ σ ≈ Time [s]
500 1000 15000100200300 u ( t =100s) σ ≈ σ ≈ σ ≈ Time [s]
600 800 1000 12000100200300 u ( t =100s) σ ≈ σ ≈ σ ≈ Time [s]
600 800 1000 12000100200300 u ( t =100s) σ ≈ σ ≈ σ ≈ ODE ODEStochastic production Stochastic production k + c < k − A B k + c > k − I I I I I
S3 Fig.
Investigating accuracy based on accumulative signaling (without protein production and degradation). (A) Regime k + c < k − with k + c = 0 . k − ( α = 100 s − , γ = 1 s − and ζ = 7). (Left) ODE model. (Right) Stochastic proteinproduction during τ b in CM and IM. (Top) Examples of time traces. (Bottom) Histograms of number of proteins producedafter 100s with standard deviation in legend based on 1000 simulations. (B) Analogous to (A) but for regime k + c > k − with k + c = 10 k − ( α = 100 s − , γ = 1 s − and ζ = 7). CM, BM and IM are shown in blue, red and green, respectively. −100 −50 0 50 100024 x 10 δ u −100 −50 0 50 1000246 δ y −100 −50 0 50 1000123 x 10 −6 δ x −100 −50 0 50 1000123 x 10 Time [ s ] δ u −100 −50 0 50 1000123 x 10 Time [ s ] δ x −100 −50 0 50 1000246 x 10 −6 Time [ s ] δ x −100 −50 0 50 1000.511.52 x 10 u −100 −50 0 50 1001234 x 10 y −100 −50 0 50 10011.0004 x ABC x 10 CMBM
S4 Fig.
Incoherent feedforward loop: Comparison of analytical results with simulations of the stochasticdifferential equations. (A) Averages of signaling rate u (left), species y from Eq. (S42) (middle) and species x from (S41)(right) as a function of time. Analytic solutions Eqs. (S32), (S43) and (12) are shown for BM in red, while a (time averaged)time-trace from a stochastic simulation using the Euler method is shown in orange (CM is almost identical and hence is notshown). (B) Corresponding variances as a function of time for k + c > k − ( k − = 6 . × s − , k + c = 10 s − ). Analytic resultsare shown in blue for CM and in red for BM; average over time (1 s ) from numerical simulations are shown in light blue for CMand in orange for BM. (C) Corresponding variances as a function of time for k + c < k − ( k − = 6 . × s − , k + c = 10 s − ).Colors same as in (B). Remaining parameters: k + c = 10 s − , k x = 10 s − and k y = 50 s − . −20 −10 0 10 203.844.24.4 x 10 u −20 −10 0 10 2010.5510.65 y −20 −10 0 10 2011.0001 x −20 −10 0 10 2011.522.5 x 10 δ u −20 −10 0 10 2000.51 x 10 −4 δ y −20 −10 0 10 2000.511.5 x 10 −5 δ x −20 −10 0 10 200.511.52 x 10 Time [s] δ u −20 −10 0 10 2000.51 x 10 −4 Time [s] δ y −20 −10 0 10 20012 x 10 −5 Time [s] δ x BAC
CMBM
S5 Fig.
Integral feedback loop: Comparison of analytical results with simulations of the stochastic differentialequations. (A) Averages of signaling rate u (left), species y from Eq. (S60) (middle) and species x from (S59) (right) as afunction of time. Analytic solutions Eqs. (S32), (S66) and (S65) are shown for BM in red, while a (time averaged) time-tracefrom a stochastic simulation using the Euler method is shown in orange (CM is almost identical and hence is not shown).(B) Corresponding variances as a function of time for k + c > k − ( k − = 6 . × s − , k + c = 10 s − ). Analytic results areshown in blue for CM and in red for BM; numerical simulations are shown in light blue for CM and in orange for BM. (C)Corresponding variances as a function of time for k + c < k − ( k − = 6 . × s − , k + c = 10 s − ). Colors same as in (B).Remaining parameters: k + c = 10 s − , k x = 10 s − and k y = 50 s − . . . . . . . . < ( u ) > / < u > CMBM . . . . A (cid:38)(cid:48) E (cid:36)(cid:3)(cid:16)(cid:3)(cid:36)(cid:48)(cid:38) (cid:39)(cid:3)(cid:16)(cid:3)(cid:41)(cid:48)(cid:37) time r ( t ) r ( t ) r ( t ) (cid:56)(cid:81)(cid:86)(cid:92)(cid:81)(cid:70)(cid:75)(cid:85)(cid:82)(cid:81)(cid:76)(cid:93)(cid:72)(cid:71) T time (cid:54)(cid:92)(cid:81)(cid:70)(cid:75)(cid:85)(cid:82)(cid:81)(cid:76)(cid:93)(cid:72)(cid:71) (cid:37)(cid:48) (cid:56)(cid:81)(cid:86)(cid:92)(cid:81)(cid:70)(cid:75)(cid:85)(cid:82)(cid:81)(cid:76)(cid:93)(cid:72)(cid:71) (cid:54)(cid:92)(cid:81)(cid:70)(cid:75)(cid:85)(cid:82)(cid:81)(cid:76)(cid:93)(cid:72)(cid:71) time time (cid:36)(cid:48) (cid:41)(cid:48) B (cid:38)(cid:48) C D (cid:37)(cid:48) (cid:36)(cid:3)(cid:16)(cid:3)(cid:36)(cid:48) (cid:39)(cid:3)(cid:16)(cid:3)(cid:41)(cid:48)(cid:37)(cid:38)
S6 Fig.
From CM (BM) to AM (FM) for multiple receptors/ion channels. (A-D) Schematic of receptor activity intime. (A) AM emerges from N unsynchronized receptors or ion channels in CM mode. (B) N synchronized CM receptors leadto a hybrid mechanism with information encoded in the frequency of broad bursts of variable duration. (C) N unsynchronizedBM receptors provide a dense series of bursts. For large N , bursts may start overlapping, leading to variable amplitudes.(D) FM emerges from N synchronized receptors in BM mode. (E) Relative variance for a system of 8 receptors with ρN synchronized and (1 − ρ ) N unsynchronized receptors, plotted for fast dynamics in the k + c < k −−