Buffering variability in cell regulation motifs close to criticality

Bistable biological regulatory systems need to cope with stochastic noise to fine-tune their function close to bifurcation points. Here, we study stability properties of this regime in generic systems to demonstrate that cooperative interactions buffer system variability, hampering noise-induced regime shifts. Our analysis also shows that, in the considered cooperativity range, impending regime shifts can be generically detected by statistical early warning signals from distributional data. Our generic framework, based on minimal models, can be used to extract robustness and variability properties of more complex models and empirical data close to criticality.

Many biological systems self-regulate their functions through bistable circuits, which have been associated to genetic [1,2] as well as growth feedbacks [3]. In particular, positive feedback loops have long been studied in systems and synthetic biology [4][5][6]; they regulate crucial functions like enzymatic activity or gene transcriptional changes during cell fate decisions [7,8]. Autoactivating positive feedback loops, simple circuit motifs promoting bistability and fine regulation of dynamical states close to self-organised criticality, are of particular importance [9,10]. Cellular heterogeneity, i.e., random cell-to-cell variations [11], can further direct transitions [12,13] and induce regime shifts between alternative stable states of gene expression or of protein concentrations [14]. Positive feedback loops with stochastic fluctuations have been observed in a variety of system including the transcription network of E. coli [15] or in the regulation of β-galactosidase [16], which results from a sudden transition from low ("off") to high ("on") level states of the lac operon at a critical point of an inducer concentration.
There are mainly two ways in which bistable systems can switch between alternative steady states [17]: transitions driven by bifurcations (which, due to loss of system resilience [18], may be anticipated by small random fluctuations) and transitions driven by large random jumps. Cells and other biological systems are hypothesized to live close to criticality to quickly respond to changing environmental conditions [19], but they should not respond to random environmental changes (noise) in order to maintain their evolutionary fitness. Close to criticality, the dynamical motifs have reduced resilience and the system can exhibit increasing variability in response to noise [20,21]. This is typical of nonlinear systems approaching a critical bifurcation and corresponds to augmented sensitivity to random perturbations and diverg- * daniele.proverbio@uni.lu ing response time, a phenomenon known as critical slowing down (CSD) [22,23].
Mechanisms to buffer variability while maintaining the critical state are thus necessary to finely regulate desired transitions [16] or to better cope with undesired shifts [18]. To this end, two main strategies can cooperate: moving the system state away from the bifurcation point, or deepening the basin of attraction to avoid random fluctuations pushing the system state to undesired attractors. These strategies correspond to changes in different environmental or regulatory conditions in cellular systems [24], allowing organisms to exploit different mechanisms to buffer variability close to criticality. In mono-stable systems, noise can be bound by the action of molecular compounds like microRNAs [25] as well as by temporal relays of signalling molecules [26,27]. For critical regimes in bistable processes, a key mechanism to buffer variability is identified here: the cooperative interactions tuning the activation function of positive feedback loops.
This buffering mechanism can be analysed by considering the simple and well-known adimensional model for stochastic autoactivating positive feedbacks [28,29]: This model describes the Michaelis-Menten kinetics of a transcriptional factor activator (denoted by x) [13,30], which activates its own transcriptions when bound to a responsive element (Fig. S1 [31]). System (1) arises from a model reduction of a two-variables genetic toggle switch, under the assumption of slow-fast timescale separation between the two variables [29]. Here, f (x) groups the deterministic terms, with steady statex (ẋ|x = 0), K is the basal expression rate, and c is the maximum production rate with critical value c 0 marking a saddlenode bifurcation (Fig. 1). The dissociation constant in the denominator of the Hill function was normalised to 1 without loss of generality [32]. The noise term η(t) accounts for intrinsic stochasticity of biological processes [33]. We consider additive Gaussian white noise to approximate the fast degrees of freedom associated to a mean field regime [21,34], with the statistical properties The Hill coefficient n, which describes the nonlinear cooperative binding mechanisms, is usually interpreted as the number of transcription factors that cooperatively promote transcription [28]. The smallest value inducing bistability in the circuit is n = 2, while n → ∞ yields the logic approximation for the activating Hill function, where Θ(·) is the Heaviside step function, making system (1) a discrete switch without bistability nor CSD. This Letter investigates the dependence of resilience properties [24] on the cooperativity index n. This way, we test how cells can keep production rates c close to their critical values and nonetheless increase resilience and buffer variability using other regulation mechanisms. As shown in Fig. 1, increasing n from 2 yields different bifurcation diagrams, where critical points shift to the left and the distance between the upper stable manifold and the unstable manifold decreases when the system gets close to criticality. We focus on systems residing on the upper branch (to be consistent with the mean field assumption) and moving left towards the saddle-node bifurcation point. This way, we investigate how biological circuits can buffer variability close to critical states by exploiting dynamical mechanisms. We also assess the parameter range where CSD-based early warning signals correctly indicate impending regime shifts.
To characterise the system stability properties, we analyse the stationary potentials and probability density functions (PDF) depending on n, in analogy to previous works [35,36]. Consider the forward Fokker-Plank where f (x) lumps the deterministic terms of Eq. (1). The stationary solution P s (x) takes the form [37] where φ(x) describes the adjoint stochastic potential whose depth is related to system resilience, i.e., its ability to recover after a perturbation. N c is a normalization constant such that Ω P s (x) = 1 (Ω is the domain). Fig. 2 shows the dependency of φ(x) and P s (x) on n when the system is either in an "off" state far away from criticality ( Fig. 2a,b), close to the criticality (Fig. 2c,d) or beyond it, where the "on" state is favoured (Fig. 2e,f). Increasing the cooperativity index n does not alter the underlying bistability, but modifies the depth of the potential and increases the separation of alternative states (Fig. S2 [31]). For the "off" and "on" states, the corresponding equilibria exhibit significantly deep attractor basins with only minor dependence on n, as also indicated in the bifurcation diagrams ( Fig. 1) and by P s (x). Close to critical points, the picture changes. The potential φ(x) displays two commensurable wells, which are more evident and symmetric for larger n, suggesting that both states become equally occupied in noisy environments. For increasing n, P s (x) displays sharper peak separation between the bistable states: the system diffuses less to intermediate states and is more constrained around single equilibrium values, as anticipated due the steeper potential barriers in φ(x). Random deviations are thus suppressed faster and transitions from one state to another are sharper and thus more robust against noise.
We now focus on how n influences variability measures, like variance and autocorrelation, close to criticality. Obtaining globally analytic expressions is challenging, in particular for high values of n. Hence, we focus on a local analysis close to the bifurcation points and employ a geometrical methodology. To derive generic results for critical manifolds, we use their local topological equivalence to bifurcation normal forms [38]. The normal forms associated to dynamical systems are simplified minimalorder forms to which all systems exhibiting a certain type of bifurcation are, around the equilibrium, topologically equivalent [39]. Supplementary Material [31] provides a brief background to normal forms and terminology used in this work. For saddle-node bifurcations like in Fig. 1, the associated normal form is [39] with two equilibrium manifoldsx 1,2 = ± √ p, one stable (+) and the other unstable (−). Note that the normal form corresponds to a parabola. To study the behaviour of stochastic solutions near the stable manifold, consider the evolution of its first-order perturbation, y = δx|x 1 exposed to the same additive white noise η(t) as in Eq. (1) [40]. Since k = 2 √ p is the distance of the control parameter value from its critical value p 0 = 0, note that k is proportional to c − c 0 from the original system, following normal form properties. The corresponding Langevin equation accounting for mean field fluctuations around the stable equilibrium is then given bẏ Eq. (7) is a typical Ornstein-Uhlenbeck (OU) process with exact solutions for statistical moments [37].
To connect the quantitative effects of n with the more qualitative topological form Eq. (6), recall that n widens or narrows the local parabolic shape of the original bifurcation diagram for Eq. (1) (Fig. 1). Eq. (6) thus needs to be augmented with a term ρ to modify the focal width of its parabolic stable manifold, which corresponds to the width of the parabola at the focal point. This leads tȯ In this formulation, ρ corresponds to the focal width of the normal form. Supplementary Material [31] contains analytical derivations for the approximation of system (1) 0.06 to the normal form (8), and its relationships with the geometrical results. Propagating ρ into Eq. (7) adds a tuning term to the bifurcation parameter, k → √ ρk. Hence, the corresponding OU process for a semi-quantitative saddlenode normal form iṡ Among its statistical moments and power spectral properties, we are primarily interested in quasi-steadystate variance (Var) and lag-1 autocorrelation (AC1), measures of system variability close to criticality. They have been proposed as proxies for system resilience and early warning signals (EWS) of impending bifurcation points [18,41]. Based on our mapping to the OU process (Eq. (9)), the analytical solutions for Var and AC1 take the form [37]: Eqs. (10) are generic for noisy saddle-node bifurcations. Fig. 3a,b shows them as functions of ρ and k. To connect with the original autoactivating feedback system, we estimate the focal width of the bifurcation diagrams for each n by fitting a parabolic form c = αx 2 + βx + γ to the data points of each bifurcation diagram in the vicinity of the saddle point. Using Matlab Curve Fitting toolbox also provides uncertainties over θ = [α, β, γ], resulting from small deviations from a perfect parabolic shape. By definition, the fitted focal width is where (x F , c F ) are the coordinates of the parabolic focus.
To get a reasonable estimate of the corresponding uncertainties, the associated standard deviation is derived from the fitted parameter uncertainties std(θ i ) using a first-order approximated propagation method [42]: The relationship between F W and n is plotted in Fig. 3c, with the corresponding std(F W ). The pattern decreases quadratically, thereby marking a rapid decrease followed by almost plateauing. Hence, a bounded and relatively small cooperativity index is in principle sufficient to effectively buffer variability close to criticality.
The estimatedρ values from fitted focal widths, for n = 2 to n = 8, are obtained aŝ where ξ is a tuning parameter proportional to the Hill function (Supplementary Material [31]). Meanρ values are marked in Fig. 3a,b with solid vertical lines. Consistently with the trend observed in Fig. 3c, the mean values spread as n increases (from left to right). Low n values yield higher sensitivity to noise, as both Var and AC1 show substantially higher values for small cooperativity indices n, even when k is large (i.e., further away from the critical point, but still within the bistable region, cf. Fig. 1). Thus, values of ρ can belong to two regions: one, where the values for both metrics are high for all k (left side of Fig. 3a,b), or another one where both metrics maintain low values for most k and increase rapidly close to criticality (right side of Fig. 3a,b). The regionρ → ∞ corresponds to the logic approximation (2) with n → ∞, where Var and AC1 also change abruptly in a step-wise manner. The ultra-sensitive regionρ → 0 is spanned by increasing dissociation constants (Fig. S3 [31]), and potentially by changing other parameters, here not explicitly considered, or by different activation functions describing, for example, wild-type vs mutant organisms [33]. Other pathways like growth feedbacks [3] will likely correspond to additional regions in the parameter space. These investigations are left to future studies. We finally investigate the performance of EWS against impending bifurcation points. The motivation is the following: consider complex systems lacking validated mechanistic models; in our case, this would translate to a scenario where n-or even the precise activation function-of an eukaryotic cell is poorly identifiable [33]. This consideration leads to questioning if we can identify statistical signals, computed on empirical data, that provide reliable information about the system's loss of resilience. Increasing trends of Var and AC1 have been widely suggested to work as EWS [18,41] but their robustness remains elusive. To study how generic they are in the identified parameter range and to account for mean trends and uncertainties, we numerically integrate the original stochastic system (1) using the Euler-Maruyama method. To mimic cell populations slowly evolving close to equilibrium, we sample 10 4 time points over 200 repeated experiments in dependence of c. This leads to a distribution of statistical indicators (e.g., see Fig. 4a, inset).
To distinguish between bifurcation-induced transitions, anticipated by loss of resilience, and noise-induced transitions, we measure the scale between the distance to the bifurcation point and the noise level by the Kramers escape rate τ = 2π( [37]. For any saddle-node bifurcation manifold (8) equipped with additive noise, Values lying at the exponential boundary of ( Fig. 3d) provide comparable ranges of control parameters and noise levels for all simulations with different n. They distinguish two regimes, one where few noise-induced transitions might occur (τ 2, Fig. 4a,b) and another regime primarily determined by bifurcationdriven resilience loss (τ 2, Fig. 4c,d). For the considered σ = 0.02, the system is very close to critical points. Fig. 4a-d display average values for Var and AC1 from numerical simulations. When the dynamics is mostly characterised by the bifurcation (Fig. 4c,d), both measures display patterns consistent with those predicted in Fig. 3a,b and increasing n better buffers variability. When the noise level becomes comparable to the potential depths, Vars for different n become very close to each other due to the more prominent role of noise-induced uncertainties (Fig. 4a). By contrast, AC1s (Fig. 4b) remain separated due to their lower sensitivity to noise (cf. Eq. (10)), but with less marked-and, therefore, harder-todetect-trends, similarly to those observed in real-world data [43].
For online applications (i.e., as new data come in and without future knowledge of the system evolution), it is necessary to quantify whether an observed increasing trend is statistically significant, assessing whether it corresponds to a EWS or some spurious fluctuation [44]. To do so, we look for significant p-values between the computed distributions close to criticality and those far from the bifurcation ("reference") (inset in Fig. 4a).  Fig. 4a,b) and farther from the bifurcation in Fig. 4f (c − c 0 corresponding to the parameters in Fig. 4c,d). The p-values cross their significance levels (either 0.1 or 0.05 [45], dashed lines in Fig. 4e,f) before the bifurcation point. This assesses that significant increasing trends of statistical indicators can be detected prior to the transitions, thus constituting reliable early warning signals. This analysis thus certifies the potential use of proposed EWS to detect approaching bifurcation points in biological motifs, providing a quantification of how much in advance the EWS become significant depending on the reference and on the p-values threshold.
Overall, our study characterised fundamental dynamical mechanisms to buffer systems' variability in critical regimes. We determined parameter ranges, corresponding to plausible cooperativity values for the positive feedback loop motif, where both variance and autocorrelation display low relative sensitivity to additive noise. In other ranges, however, the system poorly buffers its variability. Investigating whether these ranges could correspond to other dynamical mechanisms is demanded for future studies. Moreover, state-dependent noise can be further incorporated in model (1) to make it closer to biological reality [33]. Although it was not explicitly considered in this paper, primarily to focus on the effects of a single parameter on the system's stability properties, in Supplementary Material [31] (also including Refs [46][47][48]) we investigate its influence on the variability metrics considered above. Notably, it does not alter trends of AC(1) but affect those of Var (Fig. S4 [31]), as expected from dependencies on σ in Eq. (10). This is relevant in further buffering fluctuation amplitudes near criticality and calls for caution when using Var as an EWS indicator in systems strongly characterised by such type of noise. A deep investigation of interplays between noise types and resilience properties is left for future works. These may unravel alternative ways by which cells regulate their states or support the hypothesis of a self-organised finetuning in "safe" parameter spaces. Overall, our analysis contributes with quantitative insights to analytical and experimental studies of bistable systems' resilience and connects general and system-specific predictions [24].
We also assessed the sensitivity of proposed EWS to an additional regulation mechanism and suggested how to on-line quantify significant increasing trends from distributional data. We showed that extra parameters in the dynamical system do not alter the warning capabilities of indicators associated with CSD. In the considered parameters' range, they are sufficiently generic to detect resilience loss. As several indicators have been developed upon, to detect cell-fate decisions [20] and to possibly anticipate undesired shifts to, e.g., cancerous states [49,50], our results constitute an important step to interpret and apply them correctly. However, their use should be treated carefully if other quantitative mechanisms between noise and bifurcations might be at play, like certain types of state-dependent noise, possibly shadowing theoretical trends. Following our methodology, future studies might inquire other indicators and their behaviour under changing n and additional conditions. Our framework can also be easily extended to inquire the performance of buffers and EWS in different dynamical models and experimental setups.