Abstract
Background
In cell signaling terminology, adaptation refers to a system's capability of returning to its equilibrium upon a transient response. To achieve this, a network has to be both sensitive and precise. Namely, the system must display a significant output response upon stimulation, and later on return to prestimulation levels. If the system settles at the exact same equilibrium, adaptation is said to be 'perfect'. Examples of adaptation mechanisms include temperature regulation, calcium regulation and bacterial chemotaxis.
Results
We present models of the simplest adaptation architecture, a twostate protein system, in a stochastic setting. Furthermore, we consider differences between individual and collective adaptive behavior, and show how our system displays foldchange detection properties. Our analysis and simulations highlight why adaptation needs to be understood in terms of probability, and not in strict numbers of molecules. Most importantly, selection of appropriate parameters in this simple linear setting may yield populations of cells displaying adaptation, while single cells do not.
Conclusions
Single cell behavior cannot be inferred from population measurements and, sometimes, collective behavior cannot be determined from the individuals. By consequence, adaptation can many times be considered a purely emergent property of the collective system. This is a clear example where biological ergodicity cannot be assumed, just as is also the case when cell replication rates are not homogeneous, or depend on the cell state. Our analysis shows, for the first time, how ergodicity cannot be taken for granted in simple linear examples either. The latter holds even when cells are considered isolated and devoid of replication capabilities (cellcycle arrested). We also show how a simple linear adaptation scheme displays foldchange detection properties, and how rupture of ergodicity prevails in scenarios where transitions between protein states are mediated by other molecular species in the system, such as phosphatases and kinases.
Background
Chemical reactions inside cells have long been correctly described as both discrete and stochastic [13], often entailing acute spatial patterns or dependencies [46]. Despite the intrinsic uncertainty in the occurrence of these chemical events, and basically against all odds, cells prevail as efficient decision makers. Not only are their fate decisions influenced by stochastic events and embedded within widely fluctuating environments, but they are stochastic themselves [7], the underlying mechanisms of which remain widely unknown.
So, one cannot help but wonder: how do cells process widely varying information from their environment, control their own chemical 'noise', and still manage to produce appropriate responses? The key to this question lies in signal transduction pathways, a series of interconnected chemical events that lead to highly specific cell responses. One such mechanism is adaptation, a common term used to represent sets of chemical reactions that generate a transient response in the presence of a sustained stimulus [8]. These transient responses have been shown to affect gene expression and regulatory processes, where the cell decision is determined by the strength and duration of the input signal [9].
Adaptive behavior can result from three basic signaling motifs: integral control, negative feedback, and feedforward regulation [8]. The first is an abstraction of an engineering principle, where regulation is achieved by integrating the differences between a desired response and the state of the system. A cellular system may proceed in a similar fashion, by comparing 'actual' to 'desired' conditions, as has been found to be the case in bacterial chemotaxis [1013] or calcium homeostasis [14].
Integral control can be achieved through appropriate combination of negative feedback loops, the latter of which are ubiquitous elements of signaling pathways, allowing for myriads of types of physiological homeostasis. In a selfregulating gene, a transcriptional repressor negatively regulates its own expression and, within certain network architectures and ranges of feedback strength, noise can be effectively reduced. In this sense, negative feedback allows a system to respond by decreasing the magnitude of any input perturbation, generally resulting in stabilization of the input signal. However, while the latter is true in a deterministic setting, several types of nonclassic behavior can be observed once considering discrete signals and stochasticity [15].
In contrast, feedforward architectures let the system respond to known cues (input signals) in a predetermined way, independently of the system's response. This is the essential difference from feedback mechanisms, where the output influences ('feeds back') the system to create a new response. For feedforward to produce adaptation, two signaldependent pathways must affect a third component, in opposite ways, otherwise known as 'incoherent' feedforward loops [16].
Several exhaustive studies have shown that negative feedback regulation rarely yields perfect adaptation, whereas integral control and feedforward regularly do so [8,10,11,17]. Nevertheless, it should be noted that negative feedback can produce adaptation states close to 'perfect', and basically indistinguishable in terms of biological functions [8]. By perfect adaptation it is generally understood that the system will return to the exact state where it was before the input signal was introduced, provided the system was already in equilibrium.
Furthermore, in order to consider a system adaptive, certain eligibility criteria in terms of amplitude and duration of the system response have to be met. It should be noted that no homogeneous criteria exist in the literature, and comparison between different adaptation models can become a daunting task. Quite generally, though, amplitude has been assessed in terms of sensitivity and precision, namely, the difference between maximal response and prestimulation values, and the difference between equilibrium values before and after stimulation, respectively [18,19].
Recently, some types of adaptive systems (such as the incoherent feedforward loop) have been shown to display foldchange detection (FCD) properties. Namely, that the system generates a response to foldchanges in the input signal, rather than absolute levels [20,21]. The latter is related to Weber's law, which describes the relationship between a stimulus and its perceived intensity, a widely used concept in perception studies.
In this respect, some experimental studies have shown how important transduction mechanisms (such as ERK2 translocation [22] or Wnt signaling [23]) display robust foldchange responses. From these studies, several hypotheses have already arisen, such as whether cells detect and process information in relative rather than absolute terms, or whether foldchange detection facilitates the production of adjustable noise filters. Proving such hypotheses would greatly aid our understanding of cell signaling pathways, as FCD could rescale meaningful signal changes with respect to the background noise.
With all these points in mind, and in response to some of the open questions posed in [18], we study the effects of stochasticity in a minimalistic adaptation architecture, a 'twostate protein' scheme [24,25]. For such, we wanted to analyze how stochastic profiles in a singlecell system propagate to population behavior, and what this actually entails in terms of system predictability. Surprisingly, our preliminary simulations highlighted how single cell and population behavior can be completely different, adaptation largely being an emergent property of a large ensemble. This led us to analyze adaptation in an exact stochastic setting, and understand why one should think of adaptation processes in probability space, rather than in numbers of molecules.
Until now, no one had noticed how ergodicity breaks down in simple linear scenarios devoid of cell growth and replication properties. Hence, our results provide key novel insights that need to be considered in any future study of adaptation, as well as any study where biological ergodicity is readily assumed. An example of the latter is linear and nonlinear signaling pathway studies.
Additionally, we also respond to some of the open questions in [21] and show how the simple linear 'twostate protein' scheme in a stochastic setting displays foldchange detection properties, both for consecutive stimulation inputs and separate foldstimulations. This is the first study of FCD under stochasticity, the importance of which extrapolates to any cell signaling study.
Lastly, we discuss how extensions of the 'twostate protein' scheme (by considering discrete mediators, e.g. kinases and phosphatases) retain many of the properties observed in the purely linear system, including rupture of ergodicity.
Results
In what follows we shall consider the simplest adaptive architecture, previously described in the literature as a single 'twostate protein' scheme [24,25]. The idea behind this model is to consider a protein in an unmodified and modified state, denoted as P and P_{m}, respectively. The switch between the two states of the protein has basal rates k_{f }and k_{r}, and is additionally influenced by an input signal A with rate k_{a}. Moreover, the total protein concentration is allowed to vary in time, and neither the synthesis of the unmodified protein (k_{s}), nor the distinct degradation of the two states of the protein (d_{P }and ), depend directly on the input signal (Figure 1A).
Figure 1. Twostate protein models. (a) Basic linear scheme, and (b) extension to account for kinases and phosphatases mediating changes between states.
The analysis in [25] highlights several key issues. First, the steady state values of the protein states are proportional to the input signal, one being inversely proportional, the other being directly proportional. Second, when either protein degradation rate is close to zero (or much smaller than the other), the steady state concentration of the other protein is (nearly) independent of the signal. Such independence hints at the potential adaptation to changes in the input signal, a feature that can be readily observed when considering several parameter combinations. Furthermore, the choice of parameters will determine the sensitivity with respect to repeated stimulation. Namely, whether the system responds and adapts to consecutive input signals, or not.
The key questions here are: (1) which of the above mentioned properties hold when discreteness and stochasticity are taken into account, (2) do these properties only hold at the singlecell level, or extend to multicell schemes implying a kind of biological ergodicity, and lastly (3) does this system display foldchange detection properties?
Adaptation in a stochastic setting
The chemical master equation (CME) describes the time evolution of the probability P(X, t), for having x = [x_{1},...x_{N}] molecules at time t in a system with R elementary reactions, N molecular species, and volume Ω. Namely, it describes the evolution of all possible states of a chemical system in probability space.
If we consider the 'twostate protein' scheme (Figure 1A) in a stochastic setting, its time evolution will be described by the following CME:
where the index n (m) in P_{n.m }denotes the number of proteins in state P (P_{m}). Additionally, we prescribe an input signal following the step function A= S·F_{i}, where S is in units of μM
the factor F_{i }is the i^{th }signal scale parameter of F = (0.01,0.1,1,10,100), and the macroscopic reaction rate constants are set as follows: k_{s }= 0.01, k_{a }= 1, k_{f }= 1, k_{r }= 10, d_{P }= 0.01, . As is usual in a stochastic analysis, concentrations and 0^{th}/2^{nd }order reaction rates have to be scaled by the factor V·A_{v}, where V is the volume and A_{v }denotes Avogadro's constant (approximately 6.02214179 × 10^{23 }molecules^{1}). Here, we consider a volume of 1 femtoliter and selected initial values close to equilibrium. Namely P = 0.1 and P_{m }= 0.01 μM, corresponding to 60 and 6 molecules, respectively.
Now, stochastic processes can be studied by trajectory based approaches or by obtaining their underlying probability distribution function (PDF), which tracks how the probability of having specific numbers of molecules in the system changes over time. This is a daunting  many times unfeasible  task, given the combinatorial explosion of the number of coupled differential equations to consider, corresponding to increasing numbers of possible states of the system.
In fact, whenever a system is solely composed of 0^{th}/1^{st }order reactions, exact analytical PDF solutions can be obtained [26,27]. Thus, exact analytical expressions can be derived for the first two moments, which have been shown to match the solution of the system translated to a stochastic differential equation (SDE) problem [28]. Moreover, in such linear cases the first moment of the SDE solution will converge to the ODE solution, a consequence of the linearity in the drift term.
Hence, in our case, if we compare the mean of 10,000 stochastic trajectories, the corresponding deterministic solution, and the expectation of the CME reported as the sum over equal numbers of molecules of P_{m }in the exact PDF solution, it comes as no surprise that all solutions nicely match (Figure 2). As had been previously reported, the time evolution of P_{m }can show nearperfect adaptation when , while d_{P }= 0 yields perfect adaptation. However, and quite intriguingly, if we focus on individual SSA trajectories, no adaptive behavior can be inferred from single trajectories for certain values of the input signal (e.g. Figure 3C, corresponding to F_{i }= 1). The key behind this issue lies in what the solution of the CME really entails: a set of time dependent values in probability space. In this sense, computing the expected value masks how often discrepancies from this mean can happen, as well as their potential magnitude, possibly leading to entirely different dynamics. Hence, a more appropriate description of the system lies in describing the evolution of the probability, and not the expectation, in time.
Figure 2. Solutions match when using different modeling regimes (F_{i }= 1). Mean of 10,000 simulations using (A) d_{P }= 0, yielding perfect adaptation, and (B) d_{P }= 0.01, corresponding to nearperfect adaptation. Trajectories portray SSA simulations (blue), the first moment of exact CME solution using finite state projections (red), and solution of corresponding ODE system (yellow). Input signal 'stimulation' intervals were spaced by 50 units of time between one another, while each subinterval was subdivided into the following time points: [50/2^{14 }50/2^{12 }50/2^{10 }50/2^{8 }50/2^{6 }50/2^{4 }50/2^{2 }25 50]. We restricted our state space to all combinations of (P, P_{m}) such that P ∈ [0,110] and P_{m }∈ [0,30], guaranteeing our obtained PDF is closer than 99.99% to the exact solution of the CME, at all considered time points.
Figure 3. Stochastic trajectories for chemical species P_{m }in the basic twostate protein model. (AE) Single and (FJ) mean over a population of 50 cells. Subfigures correspond to different signal tuning parameters: (A, F) F_{i }= 0.01, (B, G) F_{i }= 0.1, (C, H) F_{i }= 1, (D, I) F_{i }= 10, and (E, J) F_{i }= 100. Trajectories portray SSA simulations (blue) and input signal in numbers of molecules (red).
To analyze this further, let us focus on perfect adaptation systems considering F_{i }= 1, P ∈ [0,110] and P_{m }∈ [0,30], for which we will obtain exact PDF solutions of the CME by using finite state projections [29]. This implies the solution of the CME will be both exact and complete whenever these ranges cover all possible reachable states, portraying the full probability space instead of solely describing single exact trajectories.
One can notice the deterministic adaptation value in our 'twostate' protein scheme lies close to 6 molecules of P_{m}. If we now track the evolution of the probability (i.e. the exact solution of the CME) in three separate sets: 0 to 4, 5 to 6, and 7 to 30 molecules of P_{m }it can be observed there is indeed sensitivity to the input signal and relaxation to prestimulus values, albeit in a probabilistic context (Figure 4). To understand what this entails, let us first focus on the system at time t = 50 Here, the probability of all states containing 06 molecules of P_{m }will decrease once the signal is introduced (Figure 4A), as would be expected by the sudden shift of P_{m }to higher values in the deterministic setting. In contrast, as we approach t = 100, the probability of all states containing 530 molecules of P_{m }decreases (Figure 4C). However, one should notice that, at t = 100 (and any other end of each stimulation interval), the system will revolve around states with 5 to 6 molecules of P_{m }with a probability of occurrence of roughly 32% (Figure 4B). This leaves 'a lot of room' (the remaining ~68%) for the system to be located elsewhere, as can be readily observed from single cell simulations (e.g. Figure 3C). Furthermore, of these 32%, only half belong to having 6 molecules of P_{m }in the system, which is closest to the deterministic solution. Most importantly, differences in equilibrium values entail wide variations for initial conditions of the next simulation interval, yielding the observed heterogeneity in single cell behavior.
Figure 4. Evolution of the probability using finite state projections, for separate sets of numbers of molecules of P_{m}. Here, we consider the perfect adaptation system (d_{P }= 0) with F_{i }= 1, and distinguish between: (A) 0 to 4, (B) 5 to 6, and (C) 7 to 30 molecules of P_{m}.
The description above highlights a property of any stochastic description: a system cannot be guaranteed to be in any state, and probabilistic bounds are the maximal level of information. However, one should notice the signal in this example is at least 2 orders of magnitude larger than any other component of the system. Hence, the 'twostate protein' scheme provides for a counterintuitive example where local noise overrides the effect of a strong external signal, and adaptation behavior can become an emergent property only attained at the (mean) population level. Moreover, as could be expected from the deterministic solution, the probability profile also portrays loss of sensitivity with repeated stimulation.
Population behavior
Ensemble measurements (such as those obtained from flow cytometry) display cell characteristics as distributions of values calculated over large samples of cells and, many times, distributions of cells at certain stationary states are used to infer the stochastic behavior of a single cell. By doing so, a biological version of ergodicity is necessarily implied: the percentage of the cell population in a particular state is identical to the probability to find a single cell in that state [30].
Mathematically, the ergodicity of a stochastic process is guaranteed if (i) the stochastic process is a finitestate Markov process and (ii) the stochastic process converges to a stationary state for any initial condition [1]. However, such assumptions have many times been relaxed or misinterpreted and, in such cases, erroneous conclusions with respect to single cell behavior could have been derived from population data [30].
Recently, many studies in the literature have addressed this issue. Population distributions have been based on chemical Langevin equations [3,31,32], and clever solution methods have been developed by noticing a SturmLiouville operator. For instance, the authors in [31] highlight the potential bias of population measures when cell replication rates are not homogeneous, or depend on the cell state. Similarly, the authors in [30] support the latter findings and further discuss how ergodicity breaks down whenever there are chemical interactions between cells, and when single cells display 'complicated' behavior (such as bistability).
Our exact solution of the two state protein system shows how ergodicity cannot even be assumed in minimalistic linear examples. Our linear system underlies a finitestate Markov process, is expected to relax to a stationary state, and is devoid of possible effects via cell replication and chemical interactions between cells. However, for certain input signal strengths, biological ergodicity still breaks down. In a nutshell, when the signal tuning parameter is set as F_{i }< 1, no adaptive behavior is observed in either single cells or populations. In contrast, when considering F_{i }> 1, both single cells and populations exhibit adaptation. However, when F_{i }= 1, single cells need not display any apparent sensitivity to the input signal, while both sensitivity and adaptive behavior can be observed at the population level (Figure 3C and 3H).
This brings us back to the arguments presented in the stochastic adaptation section. Namely, that adaptive behavior may not be inferred from single stochastic trajectories, since relaxation states revolving around the mean value do not necessarily represent the weighted majority in terms of probability. This can be readily observed from the values at the end states, represented as a heat map in Figure 5 implying large initial condition variations for the computation of each subsequent time interval. Actually, such wide variation can be analytically expected: it has been shown that a system composed of first order reactions that is both 'open' (i.e. including protein synthesis from source, hence violating conservation of mass) and has an initial Poisson distribution will remain a Poisson distribution at any time t > 0 [27]. Our results are entirely consistent with the theory (Figure 6). Hence, an adequate sample size and any other implications of the model have to be inferred from the underlying properties of this distribution (or else be thought as independent Binomial trials, the limit of which is the Poisson distribution).
Figure 5. Heat map of CME solution in PDF form, perfect adaptation case (d_{P }= 0). Input signal 'stimulation' intervals were spaced by 50 units of time between one another, while each subinterval was subdivided into the following time points: [50/2^{14 }50/2^{12 }50/2^{10 }50/2^{8 }50/2^{6 }50/2^{4 }50/2^{2 }25 50]. We restricted our state space to all combinations of (P, P_{m}) such that P ∈ [0,110] and P_{m }∈ [0,30], guaranteeing our obtained PDF is closer than 99.99% to the exact solution of the CME, at all considered time points. Probability values are colourcoded between values of 0 and 25%, where the end value corresponds to 25% or above.
Figure 6. PDF of adapted system and Poisson distribution. The exact solution of the CME is shown at the endpoint of all input signal stimulation intervals (blue), as compared with 100,000 randomly generated samples of the Poisson distribution (red), where the parameter λ is equal to the first moment of the CME at any of those particular time points.
Lastly, it should be noted there exist analytic results for 'open' and 'closed' systems of first order reactions, with arbitrary initial probability distributions. Hence, one can benefit from extracting an adequate sample size (or other properties) for arbitrary applications, provided all underlying reactions are unimolecular. In such cases, an appropriate PDF can be derived by the convolution of Poisson and Multinomial distributions, for any time t > 0 [27].
Foldchange detection properties
Recent studies have highlighted the peculiar capability of certain systems to respond to foldchanges in the input signal, rather than to absolute differences in numbers of molecules [20,21]. This is the essence of WeberFechner's law, which states that the maximal response to a change in signal is inversely proportional to the background signal or, in other words, that the ratio of the smallest increment and the background intensity of a signal is constant. This can be better understood by quoting an analogy used in [33]: while it is easy to understand whispered voices (increment signal) inside a quiet room (background signal), it is very hard to notice someone shouting in our ear during a Rock concert.
The recently coined term foldchange detection (FCD) implies both Weber's law and perfect adaptation. However, simultaneous application of Weber's law and perfect adaptation do not necessarily yield FCD [20], and sufficient conditions have been presented to obtain it. Namely, if a system can be described by and , where y corresponds to the output, u to the input signal, and x to the remaining chemical species, FCD can be achieved if the system is stable, shows perfect adaptation, and for λ > 0 it can be shown that f(λx, y, λu) = λf(x, y, u) and g(λx, y, λu) = g(x, y, u) [20].
It can be easily seen that the 'twostate protein' scheme does not satisfy these conditions, since:
Perhaps more intuitively, strict FCD properties could not be expected, since the 'twostate protein' scheme shows remarkable loss of sensitivity to repeated stimuli. Nevertheless, the criteria posed in [20] is only sufficient, and not necessary, so two questions are worth considering. First, to how many repetitions of the stimulus does FCD refer to? And second, to what extent are consecutive relaxations to a steady state expected to match?
In the analysis published in [21], parameter variations in an incoherent feedforward loop architecture are explored, and FCD is reported whenever the response to the two step stimuli is identical to within 10% in amplitude. By using this criteria in the 'twostate protein' scheme with foldchange input signal profiles A_{1 }= [2^{0 }2^{1 }... 2^{10}] and A_{2 }= 2*A_{1 }introduced at identical time points, FCD properties hold in two perspectives (Figure 7).
Figure 7. Foldchange detection properties in twostate protein model. Perfect adaptation twostate protein model with (a) A_{1 }= [2^{0 }2^{1 }2^{2 }2^{3 }2^{4 }2^{5 }2^{6 }2^{7 }2^{8 }2^{9 }2^{10}] and (b) A_{2 }= 2*A_{1}. Trajectories portray SSA simulations (blue), and solution of corresponding ODE system (yellow).
Both as an absolute ratio or by following the criteria used in [21], foldstimulation yields responses within 10% of the amplitude not only for any two consecutive inputs, but also for a number of consecutive steps. So, even when the system looses sensitivity to the input signal, repeated stimulation preserves foldchange detection properties. Moreover, all responses in profile A_{2 }are within 10% of the amplitude of those obtained with A_{1 }(Figure 8 and 9). Most importantly, these results hold in the deterministic regime, and the first moment of the stochastic system.
Figure 8. Relative difference in amplitude response with foldchanges in input signal. Perfect adaptation case with input profile A_{1 }= [2^{0 }2^{1 }2^{2 }2^{3 }2^{4 }2^{5 }2^{6 }2^{7 }2^{8 }2^{9 }2^{10}] μM and A_{2 }= 2*A_{1}. Relative differences were computed in three different ways: y^{1}_{j+1}/y^{1}_{j}, y^{1}_{j}/y^{1}_{j+1}, and y^{1}_{j}/y^{2}_{j}, where y^{1 }and y^{2 }are the ODE solutions using profiles A_{1 }and A_{2}, respectively, and index 'j' refers to the time point within interval 'j' where maximal response is attained.
Figure 9. Relative difference in amplitude response with foldchanges in input signal. Perfect adaptation case with input profile A_{1 }= [2^{0 }2^{1 }2^{2 }2^{3 }2^{4 }2^{5 }2^{6 }2^{7 }2^{8 }2^{9 }2^{10}] μM and A_{2 }= 2*A_{1}. Relative differences were computed in three different ways: (y^{1}_{j+1 } y^{1}_{j})/y^{1}_{j}, (y^{1}_{j } y^{1}_{j+1})/y^{1}_{j+1}, and (y^{1}_{j } y^{2}_{j})/y^{2}_{j} where y^{1 }and y^{2 }are the ODE solutions using profiles A_{1 }and A_{2}, respectively, and index 'j' refers to the time point within interval 'j' where maximal response is attained. This analysis is equivalent to the criterion used in [21].
Lastly, in nearperfect adaptation systems, consecutive inputs yield slightly different steady states. With increasing numbers of input stimuli, the equilibrium values converge to the perfect adaptation case (Figure 10). This was readily observable from the equilibrium analysis (cf. Methods) as the input signal term would dominate, yielding convergence of the steady state solutions.
Figure 10. Convergence of equilibrium values to the solution of a system with a degradation rate equal to zero. Here we considered input profile A_{1 }= [2^{0 }2^{1 }2^{2 }2^{3 }2^{4 }2^{5 }2^{6 }2^{7 }2^{8 }2^{9 }2^{10}] μM. Solutions denote a system with d_{P }= 0.01 (red) and d_{P }= 0 (blue).
Discussion
Our study of the 'twostate protein' in a stochastic setting lead us to simple yet illustrative examples on how biological ergodicity may be invalid, even in simple linear settings devoid of cell replication properties. A followup question could be: what would happen if the transitions between protein states were mediated by other molecular species in the system? In this case, one could consider discrete mediators (activators and deactivators, e.g. kinases and phosphatases) switching back and forth between the two protein states, as depicted in Figure 1B. Here, signal changes could refer to varying concentrations of the kinase and/or phosphatase.
Preliminary simulations show that all the properties studied for the simple twostate protein model can be achieved by this network topology. Such properties include perfect and nearperfect adaptation, as well as rupture of ergodicity, as exemplified by a varying kinase 'signal' in Figure 11. Qualitatively similar results can still be obtained by considering variations on the initial kinase and phosphatase concentrations (A(0), B(0) ∈ [0.01,0.1,1] μM) and single parameter values . The same holds for simultaneous variations of or while, interestingly, a completely different profile is obtained by simultaneous key variations of . For instance, if we set these parameter values to (10^{3}, 0.1, 0.1) a response that correlates negatively with changes in the signal is obtained (Figure 12). The latter can be explained by the high value of the ratio , effectively 'trapping' P molecules in a bound configuration with the kinase, hence decreasing the numbers of molecules of both P and P_{m }temporarily.
Figure 11. Twostate protein model with kinases and phosphatases as discrete mediators. (AE) Single and (FJ) mean over a population of 50 cells, with initial concentration D(0) = 0.1 μM. Subfigures correspond to different levels of the signal A, in this case the kinase: (A, F) F_{i }= 0.01, (B, G) F_{i }= 0.1, (C, H) F_{i }= 1, (D, I) F_{i }= 10, and (E, J) F_{i }= 100. Trajectories portray SSA simulations (blue), and input kinase signal (red). The following parameter values were used:.
Figure 12. Variations of key parameters in the twostate protein model, with kinases and phosphatases as discrete mediators. Columns correspond to simulations using (AD) parameter values of Figure 11; (EH) k_{a2}= k_{a3}= 0.01. Rows correspond to singlecell time courses of (A, E) P and (C, G) P_{m}; and the mean of a population of 50 cells of (B, F) P and (D, H) P_{m}. Trajectories portray SSA simulations (blue), and input kinase signal (red).
As can be intuitively expected, a similar effect in P_{m }can be obtained by considering key variations of parameters with a varying phosphatase 'signal' (data not shown). Nonetheless, the profile of P will differ from the previous case, as the reaction set is not symmetric (Figure 1B), while it should be kept in mind a separation of timescales in the protein degradation rates was necessary to achieve adaptation. Additionally, changes in the kinase/phosphatase signal are not equivalent to the minimalistic 'twostate protein' system, since these discrete mediators are treated as molecular species and, as such, their numbers can fluctuate in time. Moreover, changes to the signal intensity only refer to unbound signal molecules and, in contrast to the linear case, the total number of signal molecules within the system can accumulate in time. The latter can be readily observed in distinct chemical scenarios (such as Figure 12) where a considerable number of kinase/phosphatase molecules can be 'trapped' in a bound configuration, the molecules of which remain within the system irrespective of changes in the signal.
In our kinase/phosphatase mediated adaptation example, obtaining an exact PDF solution of the CME would be computationally expensive, if not unfeasible, given the explosion in the number of reachable states. In this case a finite state projection using Krylov subspaces, such as [34] would be preferable. Alternatively, if the mediators are highly concentrated, one could 'decouple' them from the rest of the molecular species, yielding sets of unimolecular reactions, which in turn allows for obtaining an analytic closed solution. We will leave such analysis for a separate publication, given the high discrepancy of parameter values so far encountered and the significant deviation from the main focus of this work. However, considerations such as this, or deriving alternative closure of moments of the CME to account for nonlinear terms, will greatly enhance future adaptation studies.
Conclusions
In this paper, we have studied the effects of stochasticity in a 'twostate protein' scheme, providing an explanation of what adaptation means and entails in a stochastic setting. Namely, that an adaptation profile can be achieved by calculating the first moment of the CME, but that the underlying probability distribution might be wide enough to prevent one from making definite quickanddirty assertions going from a single cell to the population level, or the other way around.
In this sense, adaptation can many times be considered an emergent property of the collective system, restricting modelers/experimentalists to obtain large samples of time courses in order to infer properties of the system as a whole. The characteristics of such samples (e.g. a minimum number of single cells to depict population behavior), are necessarily described by the underlying probability distribution corresponding to the solution of the chemical master equation. Here, we have presented the overlap of an exact solution (the CME solved in matrix form), an analytical solution (the convolution of Poisson and Multinomial distributions, which is reduced to solely the Poisson distribution in our case), as well as trajectorial and deterministic solutions. Our analysis highlights the source of variability in singlecell scenarios, explaining the cause for rupture in ergodicity in a simple linear reaction network. We have also provided a clear perspective on how systems analysis with varying inputs can be addressed.
Additionally, we have shown how the minimalistic 'two state protein' scheme displays foldchange detection properties in a stochastic setting. The latter refers to consecutive stimulation inputs and independent stimulations, and extends FCD properties to a system with nearperfect adaptation. Moreover, repeated stimulation preserved foldchange detection properties, despite loss of sensitivity to the input signal.
Lastly, we discussed extensions of the 'twostate protein' scheme by the consideration of discrete mediators (e.g. kinases and phosphatases). Our preliminary simulations show how such extensions retain many of the properties observed in the purely linear system, including loss of ergodicity.
To the best of our knowledge, this is the first time biological ergodicity has been shown to break down in a minimalistic linear architecture. The latter had been readily observed when cell replication rates are not homogeneous, or depend on the cell state, but never before in cellcycle arrested scenarios composed of solely linear nondelayed terms. Furthermore, this is the first time foldchange detection properties have been studied in a stochastic setting.
Methods
Brief explanation of the Finite State Projection (FSP) method
For the purposes of this paper, the models are both bounded and finite, so we restrict our notation to N dimensions. If we define a vector p ∈ ℝ^{n }such that each entry corresponds to the probability P(x;t) for each reachable state x, we can think of its time evolution as , where the matrix A = [a_{ij}] contains the propensities and a_{jj }= −Σ_{i≠j }a_{ij}, which basically means that each row of the matrix sums up to zero and the probability is conserved. Given an initial distribution p(0), the solution at time t is p(t) = exp(tA) p(0), where the matrix exponential is generally defined through its Taylor series expansion. If the reachable state space is large it may come in handy to consider a finite state projection [29], in which matrix A is replaced by A_{k}, a k × k submatrix of the true operator A, the corresponding indexed system states form the finite state projection and p(t_{f}) ≈ exp(t_{f}A_{k}) p_{k}(0) is the approximation to p(t) = exp(tA) p(0) at time t_{f}. An approximation can be gradually improved by adding reachable states up to a prespecified tolerance level.
In our case, the organization of the reachable states is shown in Figure 13 where (i, j) denotes the of molecules of (P, P_{m}), and corresponds to the enumerated state of the system j(M_{P }+ 1) + i + 1. M_{P }is the maximum number of molecules of P considered in the finite projection.
Figure 13. Scheme used in finite state projections.
Single cell vs. population measurements, from the simulation perspective
In general, there are two ways in which stochasticity can be considered. The 'single cell' type or the 'multiplecell' type (Figure 14). The first corresponds to studying possible states of the system in a single cell, by tracking molecular concentrations in time, once the system has reached 'equilibrium'. The second refers to comparison of multiple cells at one (or several) time point(s). In this study, preliminary stochastic simulations did not yield similar results for both types of analysis, the reason why we decided to report our results under both frameworks. In other applications, both frameworks yield similar simulation results as intrinsic stochastic noise is inherently a Markovian process and, by consequence, both cases can be considered mathematically equivalent. However, this is no longer the case in the twostate protein model presented here, which highlights how ergodicity can break down even in simple linear models.
Figure 14. Two different types of collecting sample points of stochastic simulations. (A) Independent runs (corresponding to different cells) where for each simulation the number of molecules of P_{m} is collected at the same time, T*. (B) Single run (single cell experiment) where the number of molecules of P_{m} is collected at equally spaced time steps, all of which lie beyond T*. T* is assumed to be a point in time beyond the average time it takes for a system to reach its steady state.
Solution methods and equilibrium values for the twostate protein model
The 'twostate' protein system (Figure 1) consists of six elementary reactions.
In a deterministic setting, the time evolution of such a system can be described as:
where A = A(t), and in this case A(t) is piecewise constant. The latter allowed us to use finite state projections in time intervals where A(t) = c, c ∈ ℛ.
It can be easily seen that, at equilibrium
and if , then P *= k_{s}/d_{P }and P_{m }* = k_{s}(k_{f }+ k_{a}A)/(d_{P}k_{r}). Similarly, if d_{P }= 0 we have and . As can be readily observed, when one of the degradation rates is equal to zero, the steady state solution of one molecular species is independent of the input signal.
In order to obtain transient solutions, we can formulate this as a matrix ODE problem. We shall consider here even more generalized problems, where b = g(t):
Once a fundamental solution φ(t) for the homogeneous problem has been obtained, we can use the method of variation of parameters to construct a particular solution:
As a stochastic differential equation (SDE) system with piecewise linear input signals, the 'twostate' protein system translates into:
and given the drift term is linear, the first moment of the SDE system will be identical to the ODE solution.
It is worth keeping in mind any system
with a fundamental solution yields a transient solution prescribed by
The fundamental solution will be
if matrix a_{1}(t) commutes for all times t_{1 }and t_{2 }such that t_{1 }≠ t_{2 }. Otherwise, methods such as the Magnus expansion (or Fer, symmetric Fer, Cayley, etc.) can be used to obtain a fundamental solution.
Authors' contributions
Both authors participated in designing the study. TML performed analysis of the CME and FSP simulations, AL performed SSA simulations. Both authors took part of the modeling, drafted, read, and approved the final manuscript.
Acknowledgements
Both authors would like to thank Stanly Steinberg, Kevin Burrage and Pablo Padilla for helpful comments and suggestions on the manuscript. The authors declare no funding sources for research and preparation of this manuscript.
References

Gardiner C: Stochastic Methods. A Handbook for the Natural and Social Sciences. 4th edition. Springer; 2009.

Gillespie D: A General Method for Numerically Simulating the Stochastic TimeEvolution of Coupled Chemical Reactions.
J Comput Phys 1976, 22:403434. Publisher Full Text

van Kampen NG: Stochastic Processes in Physics and Chemistry. Elsevier; 2007.

Andrews SS, Bray D: Stochastic simulation of chemical reactions with spatialresolution and single molecule detail.
Phys Biol 2004, 1:137151. PubMed Abstract  Publisher Full Text

Takahashi K, Arjunan SN, Tomita M: Space in systems biology of signalingpathwaystowards intracellular molecular crowding in silico.
FEBS Lett 2005, 579:17831788. PubMed Abstract  Publisher Full Text

Turner TE, Schnell S, Burrage K: Stochastic approaches for modelling in vivoreactions.
Comput Biol Chem 2004, 28:165178. PubMed Abstract  Publisher Full Text

Perkins TJ, Swain PS: Strategies for cellular decisionmaking.
Mol Syst Biol 2009, 5:326. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Behar M, Hao N, Dohlman HG, Elston TC: Mathematical and computationalanalysis of adaptation via feedback inhibition in signal transduction pathways.
Biophys J 2007, 93:806821. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hao N, Behar M, Elston TC, Dohlman HG: Systems biology analysis of G proteinand MAP kinase signaling in yeast.
Oncogene 2007, 26:32543266. PubMed Abstract  Publisher Full Text

Ni XY, Drengstig T, Ruoff P: The control of the controller: molecular mechanismsfor robust perfect adaptation and temperature compensation.
Biophys J 2009, 97:12441253. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yi TM, Huang Y, Simon MI, Doyle J: Robust perfect adaptation in bacterialchemotaxis through integral feedback control.
Proc Natl Acad Sci USA 2000, 97:46494653. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Alon U, Surette MG, Barkai N, Leibler S: Robustness in bacterial chemotaxis.
Nature 1999, 397:168171. PubMed Abstract  Publisher Full Text

Barkai N, Leibler S: Robustness in simple biochemical networks.
Nature 1997, 387:913917. PubMed Abstract  Publisher Full Text

ElSamad H, Goff JP, Khammash M: Calcium homeostasis and parturienthypocalcemia: an integral feedback perspective.
J Theor Biol 2002, 214:1729. PubMed Abstract  Publisher Full Text

MarquezLago T, Stelling J: Counterintuitive stochastic behavior of simple genecircuits with negative feedback.
Biophys J 2010, 98:19. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Levchenko A, Iglesias PA: Models of eukaryotic gradient sensing: application tochemotaxis of amoebae and neutrophils.
Biophys J 2002, 82:5063. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Muzzey D, GomezUribe CA, Mettetal JT, van Oudenaarden A: A systemslevelanalysis of perfect adaptation in yeast osmoregulation.
Cell 2009, 138:160171. PubMed Abstract  Publisher Full Text

Artyukhin AB, Wu LF, Altschuler SJ: Only two ways to achieve perfection.
Cell 2009, 138:619621. PubMed Abstract  Publisher Full Text

Ma W, Trusina A, ElSamad H, Lim WA, Tang C: Defining network topologies thatcan achieve biochemical adaptation.
Cell 2009, 138:760773. PubMed Abstract  Publisher Full Text

Shoval O, Goentoro L, Hart Y, Mayo A, Sontag E, Alon U: Foldchange detectionand scalar symmetry of sensory input fields.
Proc Natl Acad Sci USA 2010. PubMed Abstract  PubMed Central Full Text

Goentoro L, Shoval O, Kirschner MW, Alon U: The incoherent feedforward loopcan provide foldchange detection in gene regulation.
Mol Cell 2009, 36:894899. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

CohenSaidon C, Cohen AA, Sigal A, Liron Y, Alon U: Dynamics and variability ofERK2 response to EGF in individual living cells.
Mol Cell 2009, 36:885893. PubMed Abstract  Publisher Full Text

Goentoro L, Kirschner MW: Evidence that foldchange, and not absolute level, ofbetacatenin dictates Wnt signaling.
Mol Cell 2009, 36:872884. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

CsikaszNagy A, Soyer OS: Adaptive dynamics with a single twostate protein.
JR Soc Interface 2008, 5(Suppl 1):S4147. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Soyer OS, Kuwahara H, CsikaszNagy A: Regulating the total level of a signalingprotein can vary its dynamics in a range from switch like ultrasensitivity to adaptive responses.
FEBS J 2009, 276:32903298. PubMed Abstract  Publisher Full Text

Gadgil C, Lee CH, Othmer HG: A stochastic analysis of firstorder reactionnetworks.
Bull Math Biol 2005, 67:901946. PubMed Abstract  Publisher Full Text

Jahnke T, Huisinga W: Solving the chemical master equation for monomolecularreaction systems analytically.
J Math Biol 2007, 54:126. PubMed Abstract  Publisher Full Text

Khanin R, Higham DJ: Chemical Master Equation and Langevin regimes for agene transcription model.
Theoretical Computer Science 2008, 408:3140. Publisher Full Text

Munsky B, Khammash M: The finite state projection algorithm for the solution ofthe chemical master equation.
J Chem Phys 2006, 124:044104. PubMed Abstract  Publisher Full Text

Tozaki H, Kobayashi TJ, Okano H, Yamamoto R, Aihara K, Kimura H: Reconstructing the singlecelllevel behavior of a toggle switch from populationlevel measurements.
FEBS Lett 2008, 582:10671072. PubMed Abstract  Publisher Full Text

Sato K, Kaneko K: On the distribution of state values of reproducing cells.
PhysBiol 2006, 3:7482. PubMed Abstract  Publisher Full Text

Friedman N, Cai L, Xie XS: Linking stochastic dynamics to populationdistribution: an analytical framework of gene expression.
Phys Rev Lett 2006, 97:168302. PubMed Abstract  Publisher Full Text

Ferrell JE Jr: Signaling motifs and Weber's law.
Mol Cell 2009, 36:724727. PubMed Abstract  Publisher Full Text

Macnamara S, Bersani AM, Burrage K, Sidje RB: Stochastic chemical kinetics andthe total quasisteadystate assumption: application to the stochastic simulation algorithm and chemical master equation.
J Chem Phys 2008, 129:095105. PubMed Abstract  Publisher Full Text