For a simple model capturing some of the challenges of surviving in a fluctuating environment, consider a situation where some internal physiological quantities $\overrightarrow{P}=(P_1,\mathrm{\dots },P_N)$ must track fluctuating external variables $\overrightarrow{D}=(D_1,\mathrm{\dots },D_N)$. For example, the expression of a costly metabolic pathway would ideally track the availability of the relevant nutrient, or the solute concentration in the cytoplasm might track the osmolarity of the environment. In abstract terms, we describe these environmental pressures by the time-dependent $\overrightarrow{D}(t)$, and postulate that the organism fitness is determined by the average mismatch $-\sqrt{⟨{\sum }_{i=1}^N{(P_i-D_i)}^2⟩}$, a quantity we will henceforth call ‘performance’. Here and below, angular brackets denote averaging over time.

In this simple model, a given static $\overrightarrow{D}$ clearly defines a unique optimal state $\overrightarrow{P}$; the regulatory challenge is entirely due to $\overrightarrow{D}$ being a fluctuating quantity. The challenges faced by real organisms are certainly vastly more rich: even in the static case, the optimal behavior may not be unique, or even well-defined (optimal under what constraints?); and in the dynamic case, the future state of the environment could be affected by past actions of the organism. These considerations can add layers of complexity to the problem, but our minimal model is sufficient to focus on the basic issues of sensing, learning and responding to changing fluctuation statistics of external factors.

If $\overrightarrow{D}$ changes sufficiently slowly, the organism can sense it and adapt $\overrightarrow{P}$ accordingly. We, instead, are interested in the regime of rapid fluctuations. When changes in $\overrightarrow{D}$ are too rapid for the organism to match $\overrightarrow{P}$ to $\overrightarrow{D}$ exactly, it can rely on statistical structure. At the simplest level, the organism could match the mean, setting $\overrightarrow{P}\equiv ⟨\overrightarrow{D}⟩$. However, information on higher-order statistics, for example correlations between D₁ and D₂, can further inform the behavior and improve fitness.

To see this, in what follows, we will consider the minimal case of such structured fluctuations, namely a $N$-dimensional vector $\overrightarrow{D}=(D_1,\mathrm{\dots },D_N)$ undergoing a random walk in a quadratic potential (the Ornstein—Uhlenbeck process):

with mean $\overrightarrow{\overline{D}}$, fluctuation strength Γ, independent Gaussian random variables $\overrightarrow{\eta }$ with zero mean and unit variance, and the matrix $M$ defining the potential.

In this system, the relevant ‘fluctuation structure’ is determined by $M$ and Γ. In one dimension, Equation (1) gives $D$ a variance of $\mathrm{\Gamma }/M$. In two dimensions, denoting the eigenvalues of $M$ as ${\lambda }_{1,2}$, the stationary distribution of the fluctuating $\overrightarrow{D}$ is a Gaussian distribution with principal axes oriented along the eigenvectors of $M$, and standard deviations along these directions given by $\sqrt{\mathrm{\Gamma }/{\lambda }_1}$ and $\sqrt{\mathrm{\Gamma }/{\lambda }_2}$. Intuitively, we can think of the fluctuating $\overrightarrow{D}$ as filling out an ellipse (Figure 1C). Going forward, when we refer to learning fluctuation structure, we mean learning properties of $M$ and Γ.

If $M$ and Γ are known, the optimal strategy minimizing $⟨{(\overrightarrow{P}-\overrightarrow{D})}^2⟩$, where $\overrightarrow{D}(t)$ is set by Equation (1), can be computed exactly, as a function of the maximum allowed rate of change ${\parallel \dot{P}\parallel }^2$ (Liberzon, 2011). (If we do not constrain ${\parallel \dot{P}\parallel }^2$, the optimal behavior is of course $\overrightarrow{P}=\overrightarrow{D}$.) Briefly, the optimal behavior is to steer $\overrightarrow{P}$ toward the best guess of the expected future $\overrightarrow{D}$ (see Appendix 1, section 'Control theory calculation'). This best guess depends on the fluctuation structure, as illustrated by the comparison between Figure 1D and E for an isotropic and an anisotropic $M$.

However, in our problem, we will assume that $M$ and Γ do not stay constant long enough to be learned by evolution, and thus are unknown to the system. In this regime, it is not clear that the behavior of an $M$- and Γ-optimized system is relevant. Nevertheless, we will describe a regulatory architecture consisting of common regulatory elements that will adapt its responsiveness to the fluctuation structure of its input (‘learn’); for example, in the two-dimensional case, it will indeed develop the anisotropic response shown in Figure 1E. Moreover, we will find the steady-state performance of our architecture to be near-optimal, when compared to the theoretical ceiling of a system that knows $M$ and Γ perfectly.

The section above was intentionally very general. To discuss solutions available to cells, it is convenient to restrict the scope from this general formulation to a more specific metabolically-inspired case. From here onwards, let $D_i$ be the instantaneous demand in metabolite x_i (determined by external factors), and $P_i$ be the rate at which the metabolite is produced, with both defined in units of metabolite concentration per unit time. The number of components of the vector $\overrightarrow{D}$ now has the meaning of the number of metabolites, and we will denote it as $N_x$. The cell needs to match $\overrightarrow{P}$ to $\overrightarrow{D}$ (or, equivalently, maintain the homeostasis of the internal metabolite concentrations x_i).

The simplest way to solve this problem is via feedback inhibition. Consider first the case of a single metabolite $x$. If an accumulation of $x$ inhibits its own synthesis, a decreased demand will automatically translate into a decreased production. For our purposes, we will model this scenario by placing the synthesis of metabolite $x$ under the control of a regulator $a$ (e.g. a transcription factor), which is, in turn, inhibited by $x$ (Figure 2A). For simplicity, we will measure regulator activity $a$ directly in units of equivalent production of $x$. The dynamics of this system, linearized for small fluctuations of metabolite concentration $x$, can be written in the following form (see Appendix 1, section 'Simple end-product inhibition): 

Figure 2

Download asset Open asset

Here, we introduced P₀ with dimension of production (concentration per time) to render $a$ dimensionless. In Equation 2c, λ has the units of concentration × time, and setting $\lambda \equiv x_0{\tau }_a$ defines a time scale for changes in regulator activity. Assuming the dynamics of metabolite concentrations $x$ are faster than regulatory processes, and choosing the units so that $x_0=1$ and $P_0=1$, we simplify the equations to: 

We will refer to this architecture as simple end-product inhibition (SEPI). For two metabolites $\overrightarrow{x}=(x_1,x_2)$, the straightforward generalization is to have two independent copies of this circuit, with two regulators a₁, a₂ (Figure 2B). Denoting the number of regulators as $N_a$, we note that in the SEPI architecture, there are as many regulators as there are metabolites: $N_a=N_x$.

The architecture we will describe builds on this widely used regulatory motif, and relies on three added ingredients:

1. An excess of regulators: $N_a>N_x$;

2. Self-activation of regulators;

3. Nonlinear activation/repression of the regulators $a_{\mu }$ by the metabolite concentrations x_i.

Here and below, we use index μ for regulators ($\mu =1\mathrm{\dots }{N}_a$) and index $i$ for metabolites ($i=1\mathrm{\dots }{N}_x$).

These three ingredients, we claim, will be sufficient for the circuit to both learn higher order statistics and to use this information appropriately when matching the production to demand. It is important to emphasize that all three are readily accessible to cells. In fact, there are multiple ways to build regulatory circuits exhibiting the proposed behavior using common regulatory elements. To focus on the general mechanism rather than any one particular implementation, we will defer describing these example circuits until later in the text (Figure 6); here, we will consider a minimal modification of Equation (3) that contains the required ingredients:

This architecture bears a similarity to neural networks, and, as we will see, the familiar intuition about the value of extra ‘hidden nodes’ indeed holds. However, we caution the reader not to rely too heavily on this analogy. For example, here ${\sigma }_{\mu i}$ is a constant matrix describing how the activities of regulators $a_{\mu }$ control the synthesis of metabolites x_i.

For two metabolites ($N_x=2$) as in Figure 2C, each regulator is summarized by a 2-component vector ${\overrightarrow{\sigma }}_{\mu }=({\sigma }_{\mu 1},{\sigma }_{\mu 2})$; its components can be of either sign (or zero) and specify how strongly the regulator $a_{\mu }$ is activating or repressing the synthesis of metabolite x_i. For simplicity, below, we will choose these vectors to be of unit length. Then, each regulator ${\overrightarrow{\sigma }}_{\mu }$ is fully characterized by an angle in the $(x_1,x_2)$ plane, which allows for a convenient visualization of the regulatory systems (Figure 2D). The ${\sigma }_{\mu i}$ defines the regulatory logic of our system and does not change with time. The parameter $d\le 0$ allows us to tune the strength of the simple nonlinearity (Figure 2E); below we set $d=0$ (strong nonlinearity) unless explicitly stated otherwise. As we will show later, the learning behavior is also observed for more realistic functions such as the Hill function, but the simple piece-wise linear form of Equation (4) will help us relate the observed behavior to specifically nonlinearity as opposed to, for example, cooperativity (the Hill parameter tunes both simultaneously). Finally, the parameter κ reflects degradation and is assumed to be small: $\kappa \ll x_0$. Previously, for SEPI, it could be neglected, but here, it will matter due to the nonlinearity; for more details, see Appendix 1, section 'Simple end-product inhibition'. The parameters used in simulations are all listed in Appendix 1, section 'Parameters used in figures'.

Just like simple end-product inhibition in Equation (3), the modified system Equation (4) will correctly adapt production to any static demand (see Appendix 1, section 'Adaptation to static demand'). In the following, we will show that the added ingredients also enable learning the structure of fluctuating environments. For this purpose, we expose our system to demands $D(t)$ with fixed means (${\overline{D}}_i=1$) but a changing fluctuation structure.

To show that our system is able to adapt to different fluctuation structures, we probe it with changing environmental statistics, and show that it, first, learns these statistics, and, second, is able to make use of this information in its behavior.

For simplicity, we start with the 1-dimensional case (Figure 3A–F). In dimension $N_x=1$, an excess of regulators means we have both an activator $a_{+}$ and a repressor $a_{-}$ for the production of $x$ (Figure 3A). This is reminiscent of paradoxical regulation (Hart et al., 2012). We probe our system with changing environmental statistics by exposing it to a demand $D(t)$ with an increasing variance (Figure 3B,C). As a reminder, here and below, the mean demand is fixed at 1.

Figure 3

Download asset Open asset

Faced with a faster fluctuating input, our system upregulates both $a_{+}$ and $a_{-}$ while keeping $a_{+}-a_{-}$ constant ($a_{+}-a_{-}\approx \overline{D}=1$; Figure 3D). In this way, the two activity levels $a_{+}$ and $a_{-}$ encode both the mean and the variance of fluctuations. Crucially, the system makes use of the information it stores: The increased regulator activities allow future changes in $P$ to be faster. The system’s responsiveness, which we can define as $ℛ\equiv \frac{d\dot{P}}{dD}$, increases as $a_{+}+a_{-}$ (Figure 3E; see also Appendix 1, section 'Defining the system's responsiveness'). As a result, as shown in Figure 3F, our system is able to perform approximately equally well (after adaptation time) in each environment, unlike a system like simple end-product inhibition, which is unable to adapt its sensitivity. In summary, Figure 3D–F show that the simple architecture of Figure 3A can not only learn the statistics of environment fluctuations, but also ‘act upon this knowledge,’ effectively performing both computations of Figure 1B.

The idea of learning the fluctuation structure is perhaps clearer in dimension $N_x=2$, since the two demands can now be correlated with each other, and it seems intuitive that a system able to learn the typical direction of fluctuations (the angle α in Figure 3H) should be able to track the input better. Indeed, as we saw in Figure 1D–E, when environment fluctuations are anisotropic, the responsiveness of a well-adapted strategy must be anisotropic as well: the preferred direction must elicit a stronger response. Mathematically, the responsiveness $ℛ$ is now a matrix ${ℛ}_{ij}=\frac{d{\dot{P}}_i}{d{D}_j}$, and for a well-adapted system we expect its eigenvectors to align with the principal directions of $M$. In Figure 3G–L, Figure 4 and Figure 5, our discussion will focus on this two-dimensional case.

Figure 4

Download asset Open asset

Figure 5

Download asset Open asset

Figure 3G–L show the behavior of our system (Equation 4) with $N_a=5$ regulators (Figure 3G), exposed to an input structured as shown in Figure 3H, where we vary α. As represented pictorially in Figure 3I, we rotate the fluctuation structure matrix $M$ in Equation (1), keeping its eigenvalues ${\lambda }_{1,2}$ fixed with $\sqrt{{\lambda }_1/{\lambda }_2}=10$ (this fixes the ratio of major to minor semi-axis lengths).

With $N_a=5$ regulators, matching the mean value of $\overrightarrow{D}$ would leave $N_a-2=3$ degrees of freedom that can be influenced by other parameters (such as variance in each dimension and correlation between different demands). And indeed, changing environment statistics induces strong changes in the regulator state adopted by the system, with regulators better aligned to the input fluctuations reaching higher expression (Figure 3J; note the diagrams at the top, where dot size reflects the activity reached by the corresponding regulator at the end of each epoch; compare to the diagrams in Figure 3I). This activity pattern shapes the responsiveness matrix $ℛ$. Figure 3K plots the ‘learned angle’, defined as the direction of the dominant eigenvector of $ℛ$; we find that it tracks the stimulus angle. Finally, Figure 3L demonstrates that our architecture is able to make use of this learning, outperforming the SEPI system, whose responsiveness is isotropic and fixed.

In the previous section, we have shown by example (Figure 3) that the proposed regulatory architecture can learn the statistics of the environment. We now characterize systematically the conditions under which learning improves performance and compare our system to the theoretical performance ceiling. Note that unlike the general statement that learning correlations improves performance, the optimal performance ceiling is necessarily specific to a given model of the environmental fluctuations. Nevertheless, this comparison is informative.

The fluctuation structure in our model is defined by Γ and $M$. We first investigate the dependence of performance on Γ (Figure 4A), exposing our system to a two-dimensional input structured as in Figure 3H with $\sqrt{{\lambda }_1/{\lambda }_2}=10$ as before, $\alpha =\pi /4$, and a changing Γ.

Although the input is two-dimensional, changing Γ scales the overall magnitude of fluctuations, and the behavior is analogous to the simpler one-dimensional example shown in the first column of Figure 3. At $\mathrm{\Gamma }=0$ (static input), and by extension, for Γ finite but small, examining the steady state of Equation (4) shows that only $N_x=2$ out of $N_a$ regulators can be active. In this regime, our system is essentially identical to SEPI—the extra regulators, though available, are inactive—and in fact performs slightly worse. This is because at nonzero κ, the steady state of Equation (4) is slightly offset from the ideal state $⟨x_i⟩=1$. (While this effect can be corrected, it is only relevant in the parameter regime where no learning occurs, so we chose to keep Equation (4) as simple as possible; for additional discussion, see Appendix 1, section 'Performance penalty from the degradation term').

When Γ becomes sufficiently large, the first term in Equation (4) (proportional to the fluctuation size $\sqrt{\mathrm{\Gamma }}$) for one of the inactive regulators finally exceeds, on average, the degradation term. At this point, the system enters the regime where the number of active regulators exceeds $N_x$, and its performance deviates from the SEPI curve. Beyond this point, further changes to the stimulus no longer affect performance, as our system is able to adapt its responsiveness to the changing fluctuation magnitude (compare to Figure 3F). The threshold value of Γ satisfies $\sqrt{\mathrm{\Gamma }}\propto \kappa$; the proportionality coefficient of order 1 depends on the specific arrangement of regulators but can be estimated analytically (see Appendix 1, section 'The minimal Γ needed to initiate adaptation'). The theoretically predicted deviation points are indicated with arrows, and are in agreement with the simulation results. When a regulator in the system is particularly well-aligned with the dominant direction of fluctuations, the deviation occurs sooner, explaining the better performance of our system when the regulators are more numerous.

To better assess the performance of our system, we compare it to the theoretical optimum derived from control theory, which we represent with dotted lines in Figure 4A. For given $M$ and Γ, the family of optimal behaviors is parameterized by Control Input Power (CIP), defined as $\int {\parallel \dot{P}\parallel }^2𝑑t$. If $\overrightarrow{P}$ could react infinitely fast, it would track $\overrightarrow{D}$ perfectly, but increasing response speed necessarily comes at a cost (of making more sensors, or more enzymes for faster synthesis / degradation of x_i); constraining the CIP is thus a proxy for specifying the maximum tolerable cost. In order to compare our system with the optimal family of solutions, we compute $\frac{1}{T}{\int }_0^T{\parallel \dot{P}\parallel }^2𝑑t$ of our system at each Γ ($T$ is the simulation time), and compare to the performance of the optimally steered solution with a matched CIP; details of the calculation can be found in Appendix 1, section 'Control theory calculation'. Figure 4A demonstrates that the simple architecture we described not only benefits from matching its responsiveness to its input, but is in fact near-optimal when compared to any system of equivalent responsiveness.

It is important to note that for low Γ, the performance of the SEPI architecture also tracks the optimal curve. Throughout this work, our intention is not to demonstrate that SEPI is a ‘poor’ architecture. To the contrary, the surprising efficiency of SEPI has been noted before (Goyal et al., 2010; Pavlov and Ehrenberg, 2013), and Figure 4 similarly shows that at its own CIP, its performance is virtually optimal. The advantage of our learning-capable architecture derives from its ability to increase responsiveness when necessary, in the correct direction. Our simplified treatment of the SEPI architecture is not a strawman we seek to dismiss, but an example of a system that exhibits no learning.

Having investigated the effect of fluctuation variance (changing Γ), we turn to the effect of their correlation. Up to now, we subjected our system to a strongly correlated two-dimensional input with anisotropy $\sqrt{{\lambda }_1/{\lambda }_2}=10$ (illustrated, to scale, in Figure 1E). We will now consider a range of anisotropy values, down to anisotropy 1 (uncorrelated fluctuations, Figure 1D), keeping the variances of D₁ and D₂ constant, $\alpha =\pi /4$ as before, and $\mathrm{\Gamma }=0.05$.

The result is presented in Figure 4B. With $N_a=5$ or larger, our system is able to take advantage of the correlation, assuming it is strong enough to activate the learning mechanism. (In fact, its performance can reach values that exceed the theoretical ceiling achievable by any system that assumes the two dimensions of $\overrightarrow{D}$ to be independent, and thus must be exploiting the correlation in its inputs; see Appendix 1, section 'The system makes use of correlations in the input' and Appendix 1—figure 1). For $N_a=4$, the performance curve remains flat. This is because the four regulators are arranged as two independent copies of the system shown in Figure 3A (one $\{a_{+},a_{-}\}$ pair for each of the two inputs); this architecture can take advantage of the learned variance, but not the correlation. Finally, the SEPI architecture can adapt to neither variance nor correlation; its performance curve is also flat, but is lower. As expected, the advantage of our architecture manifests itself in environments with periods of large and/or strongly correlated fluctuations.

The model described above was a proof of principle, showing that simple regulatory circuits could learn the fluctuation structure of their inputs. Given the simplicity of our model, it is not to be expected that the exact dynamics of Equation (4) are replicated in real cells. However, the benefit of this simplicity is that we can now trace this behavior to its key ingredients, which we expect to be more general than the model itself: an excess of regulators, nonlinearity, and self-activation. In this section, we examine their role: first in our model (Figure 5), and then in more realistic circuits, relaxing our simplifying assumptions (Figure 6).

Figure 6

Download asset Open asset

In Figure 5, the parameter $d$ on the horizontal axis is the strength of nonlinearity (see Figure 2E), from perfectly linear at $d=-\mathrm{\infty }$, to strongly nonlinear at $d=0$. The vertical axis corresponds to an increasing number of regulators $N_a$, which we label as in Figure 2D; for completeness, we also include the simplest system with a single regulator co-activating both x₁ and x₂ (bottom row). Panel A examines the performance of our system as defined in Equation (4), that is, with self-activation included. In panel B, we remove self-activation by omitting the prefactor $a_{\mu }$ in front of the max function in Equation (4). The color scheme is chosen so that red indicates an improvement, and blue a deterioration, of the circuit performance relative to the SEPI architecture, which approximately corresponds to the point highlighted in Figure 5B. The difference between the labeled point and the SEPI architecture is that all models represented in Figure 5 include a small degradation term, which becomes important in the nonlinear regime. For the SEPI-like case, its effect on performance is negligible (see Appendix 1, section 'Performance penalty from the degradation term') . Performance is averaged over five angles α; see Appendix 1, section 'Parameters used in figures'.

Unsurprisingly, the performance of the simple SEPI-like architecture can be improved by adding extra regulators (pink region in Figure 5B): each new regulator allows the system to respond more quickly in a yet another direction of perturbation, with which it is ‘aligned’. However, such a strategy would have limited utility in a biological setting, since the marginal improvement per regulator must offset the cost of added complexity. The mechanism described here corresponds to the red area in Figure 5A. Importantly, in the presence of both nonlinearity and self-activation, even a single extra regulator ($N_a=3$) can already provide a significant benefit.

Figure 5A shows that in the context of our model, the reported behavior requires $N_a$ to exceed $N_x$, and $d$ to be sufficiently large. However, these ingredients are more general than the specific implementation in Equation (4). In our model, additional regulators were required because they supplied the slow degrees of freedom to serve as memory; such degrees of freedom could be implemented in other ways, for example, as phosphorylation or methylation (Barkai and Leibler, 1997). Similarly, while nonlinearity is essential (linear dynamics cannot couple to higher-order terms, such as fluctuation magnitude), its exact functional form may be changed while retaining the learning behavior (see Appendix 1, section 'Nonlinearity acts as a sensor of fluctuation variance'). Finally, the explicitly self-catalytic behavior of $a_{\mu }$ in our model is only one possible strategy for translating the stored memory into a faster response.

To demonstrate the generality of these ingredients, we constructed two circuits with very different architectures (Figure 6A,B), both reproducing the results of Figure 3C–F. These are not the only ways that the logic described above can be implemented; rather, our intention is to show that as long as we keep the key elements, we can relax our simplifying assumptions, such as the form of the nonlinearity and self-activation, while still retaining the ability to learn.

The first of these proposed circuits (Figure 6A) is based on a pair of allosteric enzymes with the toy nonlinearity of Figure 2E replaced by more realistic cooperative binding, and implements dynamics otherwise very similar to those shown above. In this circuit, the enzymes $E_{+}$ and $E_{-}$ can be in an active or inactive state: The active form of $E_{+}$, which we denote $E_{+}^{*}$, catalyzes the production of $x$; similarly, $E_{-}^{*}$ catalyzes degradation of $x$. In addition, the active enzymes can bind to molecules of the metabolite $x$ to control the self-catalytic activity. The total concentration of $E_{+}^{*}$, bound and unbound, then plays the role of the activating regulator $a_{+}$ from above ($a_{+}=[E_{+}^{*}]+[x{E}_{+}^{*}]$), while $E_{-}^{*}$ plays the role of the inhibitor $a_{-}$ ($a_{-}=[E_{-}^{*}]+[x{E}_{-}^{*}]$). The equations defining the dynamics are then:

Despite the extensive changes relative to Figure 3A, the system is still able to learn. Figure 6C compares its performance to a non-learning version with only the activating branch $a_{+}$, which is analogous to the single-activator SEPI architecture (compare to Figure 3F). For a detailed discussion of this more biologically plausible model, see Appendix 1, section 'A pair of allosteric enzymes'.

Our other proposed circuit (Figure 6B) differs significantly. Here, instead of seeking to match $P$ to $D$, the system maintains the homeostasis of a concentration $x$ perturbed by external factors. In this implementation, the production and degradation of $x$ are both catalyzed by a single bifunctional enzyme; the responsiveness of this circuit scales with the overall expression of the enzyme $E$, and larger fluctuations of $x$ lead to upregulation of $E$ due to the nonlinearity, as before. (For a detailed discussion, see e Appendix 1, section 'An architecture based on a bifunctional enzyme'.) Defining $A=a_{+}+a_{-}=[E]+[xE]$ as the total concentration of the enzyme $E$ in both its bound and unbound states, the bound and unbound fractions are described by Hill equations,

The dynamics of our system are:

Despite its compactness, this circuit is also able to learn (Figure 6D; compare to Figures 3F and 6C).

Interestingly, this particular logic is very similar to a small modification of the standard two-component signaling architecture (Figure 6E). In this architecture, the signal $s$ determines the concentration of the phosphorylated form $Y_P$ of the response regulator $Y$; the rapidity of the response is determined by the expression of the histidine kinase $X$, present at a much lower copy number. Although the signaling architecture of Figure 6E, at least in some parameter regimes, is known to be robust to the overall concentrations of $X$ and $Y$(Batchelor and Goulian, 2003), this robustness property applies only to the steady-state mapping from $s$ to $Y_P$, not the kinetics. Thus, much like in Figure 6B,a nonlinear activation of $X$ by $Y_P$ (known as autoregulation [Goulian, 2010] or autoamplification [Hoffer et al., 2001], and shown as a dashed arrow in Figure 6E) would endow this signaling system with self-tuned reactivity that learns the statistics of the input.

As the starting point of our regulatory architecture, we consider a basic form of end-product inhibition (SEPI). The environment is modeled by a time-dependent (externally induced) demand $D(t)$ for metabolite $x$ which is produced at a rate $P$ (controlled by the system); both $D$ and $P$ are defined in units of metabolite concentration per unit time. The depletion of the metabolite depends on its concentration $x$ and the demand $D$. Assuming first-order kinetics (or, alternatively, expanding a general function to linear order in small fluctuations of $x$) the dynamics of $x$ is:

Further, we consider the temporal evolution of the regulator activity $a$

By linearizing $h(x,a)$ around the stationary values $(x_0,a_0)$ we get

To examine this equation, we introduce the Fourier transforms $\tilde{a}(\omega )=ℱ[a(t)-a_0]$ and $\tilde{x}(\omega )=ℱ[x(t)-x_0]$ and get 

Equation (S4) shows that if the support of $\tilde{x}(\omega )$ is restricted to high frequencies, $\omega \gg {\lambda }_a$, then the degradation term ${\lambda }_a(a_0-a)$ in Equation (S3) is negligible. Including it would only add a restoring force, reducing the amplitude of fluctuations of $a$, and decreasing the performance of the system. Since we are interested in the regime of fast fluctuations of $x$, we choose to omit this term in the SEPI system. With ${\lambda }_x=1/\lambda$ we arrive at the dynamics used in the main text:

In the nonlinear system (Equation (3) of the main text), however, fast fluctuations of $x$ can cause the growth of $a$ (as discussed in the section ‘The nonlinearity as a sensor of the fluctuation variance’), thereby inducing slow frequency modes to its dynamics. Thus, in the nonlinear case, we cannot omit the degradation term.

The model considered in the main text modifies the SEPI architecture as follows:

Consider the case of a static input. We observe that if x₀ is set to 1, as in the main text, the presence of the degradation term causes the equilibrium point of these dynamics to be displaced away from $x_i=1$. Therefore, for a static input, the performance of this system—the degree of mismatch between $P_i$ and $D_i$, or, equivalently, the deviation of x_i from 1—is actually worse than the performance of the original SEPI.

While the case of a static input is irrelevant for the discussion in the main text, this slight offset leads to a performance penalty also for a fluctuating input. Indeed, time-averaging Equation (S7) shows that for any active regulator, we must have

where $f$ is the nonlinearity in our equation, $f(\gamma )=\max (d,\gamma )$. Clearly, we will in general again have $⟨x_i⟩\ne 1$; this is particularly obvious in the limit of small fluctuations when the nonlinearity is ‘not active’, such that $f(\gamma )=\gamma$.

This effect could be corrected by shifting x₀. In the interest of keeping our model as close to SEPI as possible, we chose not to do so: this penalty is only significant in the regime where no learning occurs, and is otherwise outweighed by the performance increase due to self-tuned responsiveness, with the additional benefit of simplifying the discussion in the main text. Even if optimizing x₀ could make the performance slightly closer to the optimal bound, this kind of fine-tuning seems irrelevant in a biological context.

In the main text, we use a measure for the ‘responsiveness’ of our system to changes in the demand $D(t)$. In this section it is shown in detail how this measure is defined. The central aim of the considered regulatory model is to match the time-dependent demand $\overrightarrow{D}$ with the regulated production $\overrightarrow{P}$. The temporal evolution of $\overrightarrow{P}$ is given by:

with

For a static demand $D_i=1$ the production relaxes to a constant value $P_i\approx 1$ (where we assumed small κ) and consequently ${\dot{P}}_i=0$. A deviation $\delta D_i$ from the static demand will lead to a change of the production $P_i$ - the larger ${\dot{P}}_i$, the faster the response to the changed demand. Therefore, we define the responsiveness of the production $P_i$ to the demand $D_j$ as ${ℛ}_{ij}=\frac{d{\dot{P}}_i}{d{D}_j}$. When assuming small fluctuations the explicit expression for the responsiveness is then given by: 

As an example we consider the one-dimensional system studied in Figure 3A-F in the main text for which the responsiveness is

where we used that ${\sigma }_{1+}=1$ and ${\sigma }_{1-}=-1$.

The calculation below follows the standard procedure known as Linear-Quadratic Optimal Control (LQC). For more details, we refer the reader to standard sources (e.g. Liberzon, 2011).

The problem we set is minimizing $⟨{(\overrightarrow{P}-\overrightarrow{D})}^2⟩$, where $\overrightarrow{D}$ follows the dynamics of Equation (1) in the main text,

We then wish to calculate the optimal way of steering $\overrightarrow{P}$. For simplicity, we will set the mean $\overline{D}=0$ (in the context of this abstract problem, this is equivalent to working with mean-shifted variables $\delta \overrightarrow{D}\equiv \overrightarrow{D}-\overrightarrow{\overline{D}}$ and similarly $\delta \overrightarrow{P}\equiv \overrightarrow{P}-\overrightarrow{\overline{D}}$). We can start by discretizing the above equation,

 where $\tilde{M}\equiv M\mathrm{\Delta }t$, and ξ has variance $2\mathrm{\Gamma }\mathrm{\Delta }t$. We seek to determine the optimal way to steer $P$; in other words, the function ${\varphi }_t(P_t,D_t)$ (‘control policy’) dictating how $P$ should be changed in a given time step:

We then can define our cost function, which combines a cost for the magnitude of u_t (how quickly we can change $\overrightarrow{P}$), and the difference between $\overrightarrow{P}$ and $\overrightarrow{D}$:

The $\varphi (P_t,D_t)$ describing the optimal behavior is the one that minimizes this cost. In order to solve for $\varphi (P_t,D_t)$, we follow standard control theory techniques and define the ‘cost-to-go’ function,

This function defines the smallest cost of all remaining steps; in particular, the total cost that we are trying to minimize is $V_0(0,0)$. The cost-to-go satisfies the boundary condition

and the following recursive relation:

Since our system is Gaussian, this recursive relation can be solved by adopting a quadratic ansatz:

Solving for the matrices $A_t$, $B_t$, $C_t$, and $Q_t$, gives us the following recursive relations:

Since our problem is to minimize the cost at steady state (known in control theory as an ‘infinite horizon’ problem, $N\mapsto \mathrm{\infty }$), we are interested in the fixed point of this mapping, specifically the matrices $A_{-\mathrm{\infty }},B_{-\mathrm{\infty }}$ to which this mapping converges when we start from ${\displaystyle A_N=B_N=C_N=\mathbb{𝟙}}$ and $Q_N=0$ (as defined by Equation (S19)).

Since $A_N$ is the identity matrix, all $A_t$ are proportional to the identity matrix as well: ${\displaystyle A_t={\alpha }_t\mathbb{𝟙}}$, where ${\alpha }_t=1+\frac{\rho {\alpha }_{t+1}}{\rho +{\alpha }_{t+1}}$. The fixed point of this mapping is $\alpha =\frac{1+\sqrt{1+4\rho }}{2}\ge 1$. Similarly, the fixed point of $B_t$ is ${\displaystyle B=[\mathbb{𝟙}-\frac{\rho }{\rho +\alpha }(\mathbb{𝟙}-\tilde{M}){]}^{-1}}$. Expressing this in terms of α only:

With these expressions, the optimal ‘control policy’ is defined by the value of $v$ that minimizes Equation S20. This defines the optimal way to change $\overrightarrow{P}$ for a given observed $\overrightarrow{D}$:

or, restoring the notations of the main text, including a non-zero $\overline{D}$:

This $u$ is the exact solution to the discrete version of the problem we considered. Since our simulations in this work use a discrete timestep, this is the form we use. Nevertheless, it is instructive to consider the small-$\mathrm{\Delta }t$, large-CIP limit such that $\mathrm{\Delta }t$ and $(\alpha -1)\mathrm{\Delta }t$ are both small compared to inverse eigenvalues of $M$. In this case we have, to first order in $\mathrm{\Delta }t$:

This leads to the following, and very intuitive, form of the optimal control dynamics:

In other words, at every step the change in $\overrightarrow{P}$ mirrors the average expected change in $\overrightarrow{D}$, with an extra term seeking to reduce their deviation. Note also that setting $\alpha =1$ (infinite CIP) corresponds to steering $\overrightarrow{P}$ directly to the expected value of $\overrightarrow{D}$ at the next timestep, as expected.

Figure 4B in the main text demonstrated that, as the fluctuating inputs become increasingly correlated, our architecture is able to outperform SEPI by an increasingly wide margin. The natural interpretation of this result is that the system is able to learn and exploit this correlation. Technically, however, one might note that this observation alone does not yet prove that our architecture is able to appropriately use the information it learned about specifically correlation. For example, it could be that strongly correlated inputs are somehow inducing a stronger increase in reactivity, causing the system to be generally faster, but without benefiting specifically from the correlated nature of its inputs.

Rather than focusing on excluding this specific scenario (which could be done by comparing the CIP values along the curves shown in Figure 4B), we will show that with a sufficient number of regulators, our architecture can perform better than the theoretical ceiling achievable by any strategy that assumes the inputs to be independent. This will formally prove that, at least for some parameters, our system’s improved performance must necessarily make use of the correlation of its inputs. Although the argument is somewhat academic in nature (we will prove our point using $N_a=25$ regulators), it is theoretically pleasing, and so we present it here.

Specifically, we consider a system subjected to inputs structured as shown in Figure 3H, with angle $\alpha =\pi /4$ so that the two inputs have the same variance. Appendix 1—figure 1 shows the performance of our architecture for several values of the number of regulators $N_a$, plotted as curves parameterized by the degradation rate κ. The degradation rate controls how large the $a_{\mu }$ can become: a high value of κ leads to lower average steady-state values of the regulator activities, causing the system to be less responsive to changes in $D$. Thus, κ can be used to set the CIP of the regulatory system, allowing us to plot these performance curves in the ‘performance vs. CIP’ axes traditional for control theory.

Appendix 1—figure 1

Download asset Open asset

Reassuringly, all these performance curves are located below the optimal control-theory ceiling computed for the true correlation structure of the input. However, the plot also shows the ‘best independent’ curve, defined as follows. Consider all possible matrices $M$ corresponding to uncorrelated inputs: ${\displaystyle M=\left(\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right)}$. Each such $M$ defines a family of control strategies (that would have been optimal if this $M$ were the true $M$ governing the dynamics of the input); this family is indexed by a parameter ρ as described above. A system following an (independence-assuming) strategy ${\displaystyle M=\left(\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right)}$ while faced with the actual (partially correlated) inputs will exhibit a certain performance $𝒫({\lambda }_1,{\lambda }_2,\rho )$ and a certain CIP, which we denote $\mathrm{CIP}({\lambda }_1,{\lambda }_2,\rho )$. With these notations, the ‘best independent’ curve is defined as

We note that the correlated-input CIP is different from the independent-input CIP that a given strategy would have incurred if faced by the input for which it is optimal. In particular, while the latter can be computed analytically, the former has to be evaluated in simulations. This makes the optimization procedure computationally costly; thankfully, the symmetry ensured by choosing $\alpha =\pi /4$ allows restricting the search to isotropic strategies ${\displaystyle M=\left(\begin{array}{cc}\lambda & 0\\ 0& \lambda \end{array}\right)}$, reducing the problem dimensionality from three parameters $\{{\lambda }_1,{\lambda }_2,\rho \}$ to more manageable two $\{\lambda ,\rho \}$.

The result is shown in Appendix 1—figure 1 as a dashed line. As highlighted in the inset, with enough regulators, our architecture is indeed able to outperform the theoretical ceiling of the best independence-assuming strategy. Although $N=25$ regulators is of course a regime irrelevant for biological applications, the aim of this argument was to formally prove a theoretical point, namely that the system as constructed must necessarily be making use of the correlation in the input signal, at least for some values of the parameters; by construction, the ‘best independent’ curve is a high bar to clear.

When framing the problem in the main text, we discussed it as consisting of two tasks, learning the fluctuation structure of the environment and ‘applying that knowledge’ (Figure 1B), and treated the two as conceptually separate. In particular, the calculation in section ‘S4 Control theory calculation’ is known as linear-quadratic Gaussian control (LQG) that assumes the correct model of the environment to already be known. This separation was done to simplify presentation, allowing us to encapsulate the control theory discussion, potentially less familiar to our audience, to a smaller portion of the narrative. In addition, the LQG calculation is simple, exact, and provides an upper bound on how well any other control strategy could perform.

Mechanistically, however, the two tasks are inseparable. In the language of control theory, our architecture implements an example of model predictive control: a strategy where the system response is governed by a continually updated internal model of the environment (here, the activity levels of the regulators, which encode the learned correlations between the inputs).

How could one distinguish a predictive vs. non-predictive control scheme in practice, when measuring all internal variables to decipher whether or not they encode an ‘internal model of the environment’ is infeasible? For our scheme, its ‘predictive’ ability manifests itself as follows. Imagine exposing our system to two inputs $\overrightarrow{D}(t)=⟨\overrightarrow{D}⟩+(\delta D_1(t),\delta D_2(t))$ which for a period of time are strongly correlated, with $\delta D_1(t)\approx \delta D_2(t)$. The learning process will drive the responsiveness matrix from one that was initially isotropic to one aligned with the correlated direction (in the notation of the main text, $\alpha =\pi /4$). Compare now the response of the system to an increase in D₁: $(\delta D_1,\delta D_2)=(a,0)$. The naïve (untrained) system would respond by increasing P₁ only, as would the SEPI architecture. In contrast, the ‘trained’ system, having learned the correlation between D₁ and D₂, responds by upregulating both P₁ and P₂ together. In this hypothetical example, our analysis predicts that deleting seemingly superfluous regulators would hinder or remove this ability (depending on the implementation, possibly without even affecting fitness in a static environment).

This is the behavior expected of a predictive controller: Under a model where $\delta D_1$ and $\delta D_2$ are strongly correlated, one expects this state of the environment to relax back to the diagonal $\delta D_1=\delta D_2$ in the near future. This kind of associative learning is experimentally measurable and, on an evolutionary scale, is well-known. Our scheme provides a mechanism for implementing this behavior on a physiological timescale. Another mechanism focusing on binary associations was previously described by Sorek et al., 2013.

In the main text, we argue that the nonlinearity in the dynamics of the regulator concentrations acts as a senor for the variance of the fluctuations. To see this, we consider the dynamics of one regulator that is controlling the production of one metabolite:

To simplify notation, we define $\gamma :=1-P/D$. Since the dynamics of $a$ are slow compared to $D$, the fluctuations of γ are on a faster timescale than the regulator dynamics. If the fluctuations of γ are small, the nonlinearity in the $\max$ function is ‘not activated’: $\max (d,\gamma )=\gamma$. In this case, the temporal average of $\max (d,\gamma )$ is zero. In contrast, if the fluctuations are strong enough, the nonlinearity is activated (see Appendix 1—figure 2). Then, the temporal average is positive, leading to an additional growth of $a$. Due to the resulting larger values of $a$, the response of the system becomes faster, making the match between $P$ and $D$ better and thus serving to decrease the variance of γ. As a result, the final average steady-state regulator concentration is reached if the system has decreased the variance of γ sufficiently by increasing the rapidity of its response.

Appendix 1—figure 2

Download asset Open asset

This argument makes it clear why the specific choice of nonlinearity is not particularly important. If $\dot{a}=f(x)$, then the static steady state satisfies $f(x_0)=0$. For a fast-fluctuating input this becomes

For any nonlinear $f$, as long as $f(x_0)=0$, the displacement of the original steady state will be determined by higher order statistics of the input. In particular, the rectified-linear nonlinearity in our equations can be replaced, for example, by a Hill function. Note that for the argument presented here, the eventual saturation of the response at large argument is irrelevant: the system will retain all the relevant behavior as long as the new nonlinearity is cup-convex in a sufficiently large vicinity of its static fixed point; see section ‘S9.1 A pair of allosteric enzymes’ for an explicit example.

The assumption $f(x_0)=0$ is not innocuous; in general, of course, the value of $⟨f(x)⟩$ is sensitive not only to the variance of $x$ (or other higher-order terms), but also to its mean, and building a system that is sensitive to specifically the variance requires adapting to the mean first. In our model, this is automatically accomplished by the underlying end-product inhibition architecture, which adapts the mean $P$ to mean $D=\overline{D}$, after which $x$ fluctuates around 1, no matter the value of $\overline{D}$.

Figure 4A in the main text includes arrows indicating theoretically derived threshold values of Γ above which our system (with a given ${\sigma }_{\mu i}$) will begin to adapt its timescale of response, deviating from SEPI in its behavior. Here, we show how this threshold value of Γ can be determined.

As discussed in the main text, at static input ($\mathrm{\Gamma }=0$) only $N_x$ out of $N_a$ regulators can be active. Consider the regulators that remain inactive in the static case, and imagine gradually increasing the fluctuation magnitude. Recall the equation for regulator activity dynamics:

After introducing ${\gamma }_{\mu }={\sum }_i{\sigma }_{\mu i}(1-P_i/D_i)$ we can write

If we chose $a_{\mu }$ as one of the regulators that remained inactive in the static case, we have ${\mathrm{\Delta }}_{\mu }<0$ at $\mathrm{\Gamma }=0$; as we begin increasing the fluctuation magnitude, the time-averaged ${\mathrm{\Delta }}_{\mu }$ will at first remain negative. The threshold Γ we seek is the one where the first (time-averaged) ${\mathrm{\Delta }}_{\mu }$ crosses into positive values. It is clear that if the fluctuations of ${\gamma }_{\mu }$ are so small that $\max (d,{\gamma }_{\mu })={\gamma }_{\mu }$ at all times, the system does not adapt. On the other hand, if the fluctuations are large enough and fast compared to the response of the system, they generate an effective additional growth of $a_{\mu }$. To first order, this additional growth term is proportional to the standard deviation $\sqrt{{\omega }_{\mu }}$ of ${\gamma }_{\mu }$. Therefore, we expect the fluctuations to cause a growth of $a_{\mu }$ if the additional growth term is large compared to κ, i.e. $\sqrt{\omega }⪆\text{c}\cdot \kappa$, with $c$ a constant of order 1.

The approximate value of $c$ can be determined using the following argument. With $d=0$, and assuming that ${\gamma }_{\mu }$ is, first, fluctuating on a fast timescale compared to ${\tau }_a$ and, second, is Gaussian with mean ${\overline{\gamma }}_{\mu }$ and variance ${\omega }_{\mu }$, we can average over the fluctuations in Equation (S28):

The system is in a stable steady state if either $⟨{\mathrm{\Delta }}_{\mu }⟩=0$ and $a_{\mu }\ge 0$ or $⟨{\mathrm{\Delta }}_{\mu }⟩<0$ and $a_{\mu }=0$. In the non-trivial first case we get the condition ${\overline{\gamma }}_{\mu }\le \kappa$. Approximating ${\overline{\gamma }}_{\mu }\approx 0$ one sees that the average growth rate $⟨{\mathrm{\Delta }}_{\mu }⟩$ is positive if $\sqrt{{\omega }_{\mu }}>\sqrt{2\pi }\kappa$, so that $c=\sqrt{2\pi }$. If this condition is satisfied, $a_{\mu }$ continues its growth until the separation of timescales between ${\gamma }_{\mu }$ and ${\tau }_a$ becomes invalid and ${\omega }_{\mu }$ decreases; this is the mechanism by which the system adapts to fast fluctuations.

The variance ${\omega }_{\mu }$ can be derived from the statistical properties of $D$. If the fluctuations of the demand $D$ are small it holds that ${\omega }_{\mu }\approx \delta {\widehat{D}}^T{\overrightarrow{\sigma }}_{\mu }\delta \widehat{D}$ where $\delta \widehat{D}$ is the covariance matrix of the stationary probability distribution of the fluctuations $\delta \overrightarrow{D}$ with $⟨\delta D_1^2⟩=\mathrm{\Gamma }\left(\frac{{\cos }^2\alpha }{{\lambda }_1}+\frac{{\sin }^2\alpha }{{\lambda }_2}\right)$ and $⟨\delta D_1\delta D_2⟩=\mathrm{\Gamma }\cos \alpha \sin \alpha \left(\frac{{\lambda }_1-{\lambda }_2}{{\lambda }_1{\lambda }_2}\right)$. The variance ${\omega }_{\mu }$ is then given by ${\omega }_{\mu }={\overrightarrow{\sigma }}_{\mu }^T\delta \widehat{D}{\overrightarrow{\sigma }}_{\mu }$ and the minimal value of Γ is set by the largest ${\omega }_{\mu }$ of the considered system.

In the main text, we proposed a simple model of metabolic regulation which highlighted the necessary properties for learning environment statistics, namely an excess of regulators $a_{\mu }$, self-activation, and a nonlinear regulation of $a_{\mu }$ by the metabolite concentrations x_i. To show how these properties can enable more realistic systems to learn the statistics of a fluctuating environment, here we present two biochemical implementations. The first of these implements dynamics nearly identical to those described in the main text, and the second, illustrated in Figure 5b, bears resemblance to two-component systems. As described in the main text, we do not necessarily expect either of these networks to be implemented in real biological systems ‘as is’. Instead, we use these to illustrate the diversity of systems that could use the logic described in this paper to learn statistics of their environment. For simplicity, we consider the one-dimensional case (first column of Figure 3 in the main text).

S9.1 A pair of allosteric enzymes

The circuit is shown in Appendix 1—figure 3. The enzymes $E_{+}$ and $E_{-}$ can be in an active or inactive state: The active form of $E_{+}$, which we denote $E_{+}^{*}$, catalyzes the production of $x$; similarly, $E_{-}^{*}$ catalyzes degradation of $x$. In addition, we postulate that active enzymes can bind to molecules of the metabolite $x$, which controls self-catalytic activity (see diagram). The total concentration of $E_{+}^{*}$, bound and unbound, plays the role of the activating regulator $a_{+}$ from the main text ($a_{+}=[E_{+}^{*}]+[x{E}_{+}^{*}]$), while $E_{-}^{*}$ plays the role of the inhibitor $a_{-}$ ($a_{-}=[E_{-}^{*}]+[x{E}_{-}^{*}]$).

Appendix 1—figure 3

Download asset Open asset

The same regulatory structure could be realized with transcription factor regulation, with the role of the active enzymes ($E_{+}$ and $E_{-}$) played by two transcription factors. In this version, activation/deactivation of enzymes is replaced by the simpler process of transcription factor synthesis/degradation. For concreteness, here, we focus on the enzymatic case, largely because we expect architectures like ours to be more relevant in fast-responding circuits, which tend to be non-transcriptional. However, except for the difference of timescales, the dynamics of the two versions would otherwise be identical; in particular, both implementations would ‘learn’ in the same way.

For this discussion, it is convenient to have timescales of dynamics of $a_{\mu }$ and x_i encoded as explicit parameters. Assuming first-order kinetics, the dynamics of the network can then be described by:

(S30) ${\displaystyle \{\begin{array}{rl}{\displaystyle {\tau }_x\dot{x}=}& {\displaystyle {\gamma }_{+}{a}_{+}-{\gamma }_{-}x{a}_{-}-xD(t),}\\ {\displaystyle {\tau }_a{\dot{a}}_{+}=}& {\displaystyle a_{+}\frac{c_{+}^n}{c_{+}^n+x^n}-a_{+}{\kappa }_{+},}\\ {\displaystyle {\tau }_a{\dot{a}}_{-}=}& {\displaystyle a_{-}\frac{x^m}{c_{-}^m+x^m}-a_{-}{\kappa }_{-}.}\end{array}}$

Here, we assume that the metabolite $x$ is much more abundant than the active enzymes $E_{+}^{*}$ and $E_{-}^{*}$, meaning that the relative amount of bound $x$ is very small. This allows us to neglect, in the dynamics of $x$, the source and sink terms due to binding and unbinding of $x$ to the enzymes. We also assume that this binding and unbinding occurs on a much faster timescale than all other processes.

Appendix 1—figure 4 shows an example of simulation results for these dynamics (for the full list of parameters used, see section ‘S11 Parameters used in figures’). We see that the system reacts to an increasing variance of environmental fluctuations (A) by increasing regulator activities (B). The figure also shows the behavior of a SEPI system which only consists of one $a_{+}$ regulator described by the dynamics in Equation (S30). Appendix 1—figure 4C shows that the response strength, defined as discussed in the SI section ‘Defining the system’s responsiveness,’

(S31) $ℛ=\frac{d\dot{P}}{dD}\approx {\gamma }_{+}{a}_{+}\frac{n{c}_{+}^n}{{(c_{+}^n+1)}^2}+{\gamma }_{-}{a}_{-}\frac{m{c}_{-}^m}{{(c_{-}^m+1)}^2},$

is increasing due to the changed regulator activities. Finally, Appendix 1—figure 4D compares the performance of the system Equation (S30) with the corresponding SEPI system (which, again, we define by the same equations as Equation (S30), except without the $a_{-}$ regulator). Similar to Figure 3F in the main text, the performance of the adapting system does not change as the variance of the stimulus increases, while the SEPI system becomes worse.

Appendix 1—figure 4

Download asset Open asset

S9.2 An architecture based on a bifunctional enzyme

For the reader’s convenience, we reproduce this circuit in Appendix 1—figure 5 (identical to Figure 6B in the main text). As described in the main text, for greater generality, we will here rephrase the task: instead of matching production to demand, we will think of maintaining the homeostasis of a quantity $x$ perturbed by external factors. For example, instead of being a metabolite concentration, $x$ could be osmolarity mismatch, and our circuit a hypothetical architecture for active control of osmotic pressure. In this interpretation, the enzyme $E$ might be a mechanosensor triggered by tension in the membrane or cell wall, while ‘production’ and ‘degradation’ of $x$ could be activities of opposing pumps, or regulators of glycerol export or metabolism.

Appendix 1—figure 5

Download asset Open asset

To simplify our language when describing the terms in the specific equations we simulate, we will still talk of a metabolite $x$ being produced and degraded. However, to better accommodate alternative interpretations, the regulator activities will now be defined so that $a_{+}$ and $a_{-}$ would be equal on average (for example, the activities of pumps in opposite directions should, on average, balance). This homeostasis-maintaining formulation is in contrast to Figure 3D in the main text, where regulators satisfied the constraint $⟨a_{+}-a_{-}⟩=\overline{D}=1$.

The production and degradation of $x$ are catalyzed by a bifunctional enzyme that changes its activity when bound to $x$ forming the compound $xE$. The concentration of the unbound form $E$ corresponds to the activating regulator, $a_{+}=[E]$, and increases the production $P$ of $x$, while $xE$ plays the role of the inhibiting regulator, $a_{-}=[xE]$, and promotes the degradation of $x$.

As above, we assume first-order kinetics for the production and degradation of $x$, and that the binding kinetics are fast compared to the other timescales in the problem. Defining $A=a_{+}+a_{-}=[E]+[xE]$ as the total concentration of the enzyme $E$ in both its bound and unbound states, the bound and unbound fractions are described by Hill equations:

(S32) $a_{+}=A\frac{c^m}{x^m+c^m},a_{-}=A-a_{+}=A\frac{x^m}{x^m+c^m}.$

These expressions make it explicit that a small change in the concentration $x$ induces a change in $a_{+}$ and $a_{-}$ that is proportional to their sum, $A=a_{+}+a_{-}$. In this way, the circuit shown in Appendix 1—figure 6 does include an element of self-activation (one of our ‘key ingredients’), even though no interaction in the diagram is explicitly self-catalytic.

Appendix 1—figure 6

Download asset Open asset

With these notations, the dynamics of our system are:

(S33) ${\displaystyle \{\begin{array}{rl}{\displaystyle {\tau }_x\dot{x}=}& {\displaystyle P_0+{\gamma }_{+}{a}_{+}-{\gamma }_{-}x{a}_{-}-xD(t)}\\ {\displaystyle {\tau }_A\dot{A}=}& {\displaystyle -A\kappa +f(x),}\end{array}}$

where we assumed that modifying the enzyme $E$ does not significantly affect the quantity $x$. (In the metabolic formulation, this corresponds to the assumption that $x$ is much more abundant than $E$, so that the sequestration of $x$ by $E$ has negligible effect on the free concentration of $x$; in the osmotic formulation, the act of triggering a mechanosensor has negligible instantaneous effect on pressure). In the second equation, the synthesis of the enzyme $E$ depends nonlinearly on the metabolite concentration $x$. The specific form of nonlinearity does not significantly affect the results, as long as it is sufficiently cup-convex in the vicinity of the operating point: As described in the section ‘S7 Nonlinearity acts as a sensor of fluctuation variance’, we can think of the nonlinearity $f(x)$ as a ‘sensor’ for the variance of environmental fluctuations. Whenever fluctuations in $D(t)$ increase such that the current responsiveness of the system fails to maintain the homeostasis of $x$ within previous bounds of fluctuation magnitude, the fluctuations of $x$ will lead to growth of $A$, increasing the responsiveness until it is again able to reduce the fluctuations of $x$. In our simulations we chose a Hill function with cooperativity 4 (see section ‘S11 Parameters used in figures’).

Appendix 1—figure 6 shows simulation results for this system. As in the first column of Figure 3, the variance of $D$ is increased and the response of the system to this change is monitored. We see that the regulator concentrations correspondingly increase, causing a larger response strength $|\frac{d\dot{x}}{dx}|\approx 1+\frac{2\gamma E{c}^{m}m}{{(1+c^m)}^2}$. The increase in response strength is able to compensate for most of the performance loss, which shows that the system successfully adapts its timescale of response. This is in contrast to the ‘SEPI-like’ system with a fixed value $A=1$, which cannot adapt its responsiveness and whose performance drops with every increase in fluctuation variance.

In the main text, we argue that the production $\overrightarrow{P}$ of the proposed system Equation (3) adapts to any static demand $\overrightarrow{D}$. The full dynamics of the system is

With a static demand, Equation (S34) possesses the same fixed points as the simplified dynamics:

These dynamics have a Lyapunov-function

This can be verified by considering the temporal change of $F$

with ${\mathrm{\Delta }}_{\mu }={\sum }_i{\sigma }_{\mu i}(1-P_i/D_i)-\kappa$. Thus, $F$ always increases and is obviously bound from above. For small κ, the maximum of $F$ is reached for $\overrightarrow{P}\approx \overrightarrow{D}$, showing that the system adapts to any static demand.

If not stated otherwise, we use the following parameters: $\overline{D}=1$, $\mathrm{\Gamma }=0.35$, $d=0$, $\kappa =0.015$, ${\tau }_a=3$, $\alpha ={45}^{\circ }$, ${\lambda }_1=8.75$, ${\lambda }_2=875$, $dt=1/{\lambda }_2\approx 1.14\cdot {10}^{-3}$. Since the demand $D_i$ is modeled by a stochastic process which is, in principle, unbounded, there is a non-zero probability that the demand $D_i$ becomes negative. To prevent this behavior in the simulations we set $D_i=0.01$ if $D_i<0.01$.