Information propagation in time through allosteric signaling

Naively, one expects the Shannon information communicated by a protein downstream to a signaling network, in which the protein is embedded, to grow monotonically with this protein's rate of product formation. Here we observe that this does not hold for allosterically regulated proteins often found in signaling networks. In particular, we find that the Shannon information communicated by an allosterically regulated protein is a complex function of the transition kinetics between its different allosterically regulated states and their respective product formation rates. Surprisingly, our work implies that allosteric down-regulation of protein product formation may not only be used to silence proteins as one would normally expect. Rather, down-regulation may be used to increase the Shannon information communicated by this protein in a signaling pathway.

Here, we propose a framework to address how much information is communicated by a protein (the "sender") to other proteins ("receivers") in the same signaling pathway. Naively, we may think of information communicated to receivers of a transmitted message as growing with the concentration of product generated by the sender. This naive picture is immediately complicated by allostery, a process of remote regulation of functional activity by which biological macro-molecules, such as proteins, transmit the effect of ligand binding from one site to a distal functional site [4,[10][11][12][13][14]. Several studies have been performed using information theory to analyze the role of allosteric regulation in signaling pathways [15,16]. These studies primarily explored the thermodynamics of substrate binding through allostery.
Rather, here we focus on the communication between signaling proteins through product formation events. As a result of allosteric up-or down-regulation of its functional site [10,11], a sender's product formation rate varies in time. Put differently, the Shannon information of a probability distribution, P , over the sequence of product arrivals now depends on the kinetics of transitions between the different allosteric states whose product formation rates differ. These kinetics are encoded in the sequence of product arrivals, n 1:M , indexed in discrete time from 1 to M .
In what follows, we first consider a simple master equation model capturing the essence of allostery. We compute the probability of product arrival, P (n 1 , n 2 , . . . , n M ), to the receiver (by marginalizing over allosteric states invisible to the receiver) using tools from Hidden Markov Models (HMMs) [17,18]. We then compute the Shannon information over this distribution P (n 1 , n 2 , . . . , n M ) which encodes both concentration and allosterically induced state-switching dependencies.
Allostery, first articulated over 50 years ago by Monod [10], has remained a central point of focus in biology. Since then, it has been shown that allosteric effects propagate beyond the molecular level and play a role at the level of cellular signaling [12]. Our work now suggests a broader role for allostery: a way to increase information communicated in signaling pathways even, counter-intuitively, in down-regulation. MODEL We start with a minimal model of allosteric regulation in proteins inspired by Refs. [10,11]. We consider a sender with catalytic (i.e., functional) site A which binds and unbinds substrate S to make a complex, AS, with rates k A+ and k A− respectively. This complex degrades to produce product, P , with rate d p releasing this protein's catalytic site back to its original unbound form A. Briefly, we write In addition, S can also bind an allosteric regulator site B on this same protein with association and disassociation rates k B+ and k B− respectively. The effect of BS may be allosterically propagated to site A with rate h p where, A * is the unbound catalytic site in its altered (upor down-regulated) state. In such a state, the dynamics of interaction between the catalytic site and the substrate changes as follows where, k * A+ and k * A− are, as before, association and dissociation rates and r p is the allosterically altered product formation rate. For sake of generality, we allow the catalytic site to relax with rate s p For illustrative purposes, the reactions described in Eqs. 1,2,3,4 and 5 can be straightforwardly simulated using the Gillespie algorithm [19]. We observed "bursts" of higher rates of product formation events as anticipated in the case of allosteric up-regulation, Fig. 1. Once the steady state is achieved, the waiting times between product arrivals without allostery, are approximately exponential distributed (as expected). Put differently, for the choice of parameters specified in the caption of Fig. 1, the ratio of the mean squared to the variance in product arrival is ≈ 1.00 when allosteric regulation is absent (Fig. 1b, as expected for Poisson processes) but strongly deviates from unity (≈ 17.40) otherwise (Fig. 1a). This simple result already preliminarily suggests that the information communicated in time by allostery warrants attention.
At steady state, i.e., when the bound complex concentration (AS, BS and A * S) and the product formation rate become time independent, the system of chemical reactions is simplified to where, k f = k A+ dp k A− and k * f = k A * + rp k A * − are the effective product formation rate constants for the catalytic site in states 1 and 2 respectively; and [BS] eq is the bound We observe that the model with allostery exhibit a "bursty behavior", whereas this is not seen in the model without allostery. The values of the parameters used are: kA+ = 50s −1 , kA− = 25s −1 , dp = 0.9s −1 , kB+ = 50s −1 , kB− = 25s −1 , kA * + = 75s −1 , kA * − = 25s −1 , rp = 50s −1 , sp = 70s −1 , hp = 50s −1 for (a) and 0s −1 for (b).
complex allosteric regulator concentration at steady state (see Supplementary Information). The now linear switching dynamics of the catalytic site between its two states can straightforwardly be described using transition matrices [20] P (s t+δt , s t ) = Γ t→t+δt P (s t ) where, P (s t ) is the probability vector describing the catalytic site's state at time t (state 1, first index; state 2, second index); Γ t→t+δt is the transition matrix whose elements describe the transition probabilities between states from time t to t+δt; and P (s t+δt , s t ) is a vector representing the joint catalytic site state occupancy probability at time t and t + δt. The matrix elements of Γ t→t+δt depend on rate constants and [BS] eq (see Supplementary information).
In order to calculate the information communicated by an protein, we keep track of the number of products generated per time interval, δt. Here we assume that δt is si si-1 si+1 FIG. 2. Graphical Model describing a two state system with Poissonian emission. Here the observable, ni, is the number of products produced by a catalytic site at time level i in the hidden state si. Γi−1→i represents the transition matrix of the states between time steps i − 1 to i.
small enough such that the probability of state transition is small. In this way, we can map our problem given by Eqs. 6,7 and 8 onto an HMM (details in Supplementary  Information), see Fig. 2.
The model in Fig. 2 describes the dynamics of a catalytic site existing in state s i at time step i and produces products n i as described by the Poisson process where λ si is given by (see Supplementary Information) Within this HMM framework, the probability of observing a given trajectory in the number of products in M time steps, i.e., n 1:M = {n 1 , · · · , n M }, is where, the last term, P (s i |s i−1 ) carries the dynamics of switching between states for the catalytic site according to Eq. 9, and P (s 0 ) represents the initial conditions. We call n 1:M a "product time series".
The summation over the right hand side of Eq. 12 is performed over all possible states realized at all time steps [21] (see Supplementary Information).

RESULTS AND DISCUSSION
In our model, a signal is composed of a trajectory of product formation events observed in time. As such, we define the information carried by such signals as the Shannon information [1] over the distribution of the number of products that could be observed in the time series, i.e., where the expression is summed over the entire sample space of trajectories, n i ∈ 0, 1, 2, ..., ∞ for all i ∈ [1, M ]. Since the expression for the probabilities cannot be calculated analytically, the summation also does not assume an analytic form. For this reason, to calculate the information we performed the summation over a subset of trajectories sampled by a discrete stochastic simulation of a two state HMM. That is: 1) we sample some trajectories from a discrete stochastic simulation; 2) we calculate the probabilities of such trajectories using the forward filter; 3) we calculate the information using Eq. 13 by summing over these sampled trajectories. Finally, 4) Steps 1-3 are repeated until the information converges to a pre-defined precision. Using steps 1-4 (details in Supplementary Information), we analyzed the effect of allosteric regulation on the information contained in the distribution P (n 1:M ). That is, the information encoded in product arrivals. In our model, modulation of information carried by product arrivals can be achieved through allosteric regulations by changing the average number of products generated by the catalytic site in the two states or by altering the switching probabilities between two states.
We first investigated the effect of varying the product formation rate of the first state, λ 1 by initiating the protein in the first state, whilst keeping the product formation rate of the second state as well as the switching probabilities between states fixed, Fig. 3. At first, as a sanity check, we observe that without any allosteric regulation, i.e., when the second state is never visited (P 1→2 = 0), the information increases linearly with higher values of λ 1 as one would expect from a catalytic site generating products distributed according to a Poisson distribution. Next, when λ 1 < λ 2 , i.e., in the case of up-regulation, the information communicated by the sender increases but not monotonically. That is, counter-intuitively, the information communicated by the sender does not naively grow by simply increasing λ 1 . Indeed decreasing λ 1 well below λ 2 may increase the information communicated by the sender. As expected, the effect of allostery vanishes as λ 1 = λ 2 . Equally counter-intuitively, the information communicated by the sender rises again when λ 1 > λ 2 , i.e., for the case of down-regulation, as temporal kinetics now increase the information content in the distribution of product arrival.
Second, we analyzed the variation of the information transmitted with allosteric regulation by varying the switching probability P 1→2 holding everything else fixed, Fig. 4. In the absence of regulation, i.e., λ 1 = λ 2 , the effect of allostery vanishes and information conveyed is insensitive to the value of P 1→2 . For an up-and downregulated proteins, the catalytic site conveys minimum information when restricted to state 1. On the other hand, as P 1→2 grows, with up-regulation, the catalytic site is able to access a state with higher average production and therefore conveys a larger amount of information. However, unexpectedly, with down-regulation, the catalytic site is still able to transmit a larger amount of information for a range of switching probabilities. This counter-intuitive observation illustrates once more that information is not a mere function of amount of products sent but also depends on switching kinetics between states. We observe that without any allosteric regulation, the information encoded in the product time series do not vary with switching probabilities. However, with up-regulation, i.e., when λ1(60) < λ2(100), the protein communicates more information as the probability of it switching to an up-regulated state increases. Counter-intuitively, a down-regulated protein, i.e., when λ1(100) > λ2(60), is still able to communicate a larger amount of information for a wide range of switching probabilities. For all plots, we fix P2→1 to 0.08.

CONCLUSION
Here we proposed a minimal model capturing salient allosteric features. Using this model, we quantified the information communicated by allosterically regulated proteins through time-varying product formation rates. Our choice of two state models is a matter of convenience; our formalism readily generalizes to more complex kinetic models.
More importantly, we found that an allosterically regulated protein may convey a larger amount of information as compared to a protein with no allosteric regulation by sometimes decreasing its total product output. This suggests that allostery may provide a means to control the information communicated in product time series in a way that goes beyond the energetically demanding "more product, better signal" exploitative paradigm. This is, allostery may provide "lower signal but more information". The possibility of parameter fine-tuning to communicate more information is especially relevant given allostery's key role in protein evolution [12,22,23]. It opens the possibility that Nature may fine-tune allosteric parameters (including switching rates between states as well as production rates) to adapt/evolve its signaling pathways in response to external stimuli warranting exploratory (high information/low signal) or exploitative (low information/high signal) strategies [24].
In the main text, we proposed a minimal model where a set of coupled chemical reactions describe the interaction of a catalytic site A-binding substrate S to form product P -whose activity is modulated by another binding site B through allosteric regulation.
The generation of products through this set of coupled chemical reactions are realized using Gillespie algorithm [1]. Firstly, as a proof of concept, we analyze the rate of product formation in our model for a protein with no allosteric regulation between catalytic sites A and B, i.e., when h p = 0s −1 , Fig. 1a. We observe that, as expected, without any allosteric regulation, the system of chemical reactions reduce to Michaelis-Menton kinetics [2] such that the rate of product formation increases linearly with smaller concentrations of the substrates and then saturates to a constant value at higher concentrations of substrate. On the other hand, when allosteric regulation is present, i.e., h p > 0s −1 , the production rate exhibits a non-linear relationship at lower concentrations of substrate, commonly observed in allosteric proteins [3][4][5], Fig. 1b.
As discussed in the main text, the model proposed above produces a unique "bursty" emission of products when allosteric regulation is present. The model also captures this "bursty" production arising from different microscopic models of allosteric regulation [6]. That is, namely K-type models (when allosteric regulation manifests itself by modulating the catalytic site's binding affinity with the substrate, Fig. 2) and V-type models (when allosteric regulation affects the rate by which the bound complex of catalytic site reduces to give products, Fig. 3).
Coupled chemical reactions for allosteric regulation reduce to a two state Hidden Markov model When the system reaches steady state, i.e., when the rate of product formation reaches a constant value, where [AS] eq and [A * S] eq are the steady state concentrations of the bound complexes. These concentrations are attained when the rate of forward and backward reaction for complex formation of the substrate with the catalytic site and the allosteric regulator site are equal. Also, at this state the concentrations of AS and A * S does not change with time. Therefore, at steady state, we have  Variation in rates of product formation with substrate concentration. We observe that the model without allosteric regulation, (a) recapitulates the dynamics commonly observed in proteins following Michaelis-Menton Kinetics. On the other hand, when allosteric regulation is present, (b) the rate of product formation is nonlinear at lower concentrations of substrate, a signature of allosteric regulation. The following parameter values were used: kA+ = 50s −1 , kA− = 25s −1 , dp = 5s −1 , kB+ = 50s −1 , kB− = 25s −1 , kA * + = 75s −1 , kA * − = 50s −1 , rp = 10000s −1 , sp = 50s −1 , hp = 0s −1 for (a) and hp = 100s −1 for (b). The system was simulated for 30s and rates for each substrate concentration were averaged over 1000 independent runs. as shown below Combining Eqs. 10 and 11 with Eqs. 7, 8 and 9 with the condition that substrate is present in excess, i.e.
products from catalytic site in the two states (say 1 and 2 respectively). Now, for a catalytic site in state 1, the probability of remaining in state 1, P 1→1 , is proportional to [A] eq [S] eq k f from Eq. 12. On the other hand, the probability P 1→2 is proportional to [A] eq [BS] eq h p , Eq. 14. Similarly, for a catalytic site in state 2, P 2→2 is proportional to [A * ] eq [S] eq k * f , Eq. 13 and P 2→1 is proportional to [A * ] eq s p , Eq. 14. These probabilities can be explicitly calculated using appropriate normalization conditions, i.e., using P 1→1 + P 1→2 = 1 and P 2→1 + P 2→2 = 1. These can be used to calculate the transition probability matrix for a catalytic site at steady state as:  Bursts of products observed in V-type allostery. We observe that the model with V-type allosteric regulation, (a) shows a "bursty" formation of products not observed in the model without allosteric regulation, (b).

COMPUTING THE PROBABILITY OF OBSERVING A TRAJECTORY OF PRODUCT FORMATION EVENTS WITH ALLOSTERIC REGULATION
As shown above, the minimal model for allosteric regulation can be reduced to a two state model. In this two state model, we observe product formation events in time. In order to further simplify the model, we move from a continuous time regime to a discrete time regime where the observable is the number of product formation events detected in a time level.
Upon discretizing time, the time interval, δt is selected such that the probability of observing the product formation events from a catalytic site in a given state is much larger than the probability of switching the state of a catalytic site. This ensures a separation of time scales between the product formation events and the state switch-ing events as shown below (16) and, This enables us to represent the model as a Hidden Markov Model where a catalytic site stochastically switches states quickly between two time levels which are relatively larger. This scheme can be shown using the transition matrix, Γ, described earlier, where the offdiagonal elements representing the probabilities of the catalytic site switching its state will be much smaller than the diagonal elements. Therefore, the dynamics of switching from state s i at time level i to state s i+1 at time level i + 1 can be described as: Here, P (s i ) represents the probability distribution of the state occupied by the catalytic site at a time level i. Moreover on account of our steady state assumptions, the elements of the transition matrix, Γ remain constant as described in Eq. 15, therefore, Γ i+1→i = Γ: ∀i ∈ 1, 2, ..., M . The probability of observing n i bursts of products by a catalytic site in state s i in time level i is where λ si is the average number of products formed by a catalytic site in state s i which depends on the steady state concentration of the bound complex in that state s i as well as the rate parameters as follows The probability of producing a series of n 1:M products while the catalytic site spends M time levels in states where, each s i ∈ (1, 2), is the state of the catalytic site at time level i, and P (s 0 ) represents the initial condition of the catalytic site. Therefore, for a catalytic site initialized in state 1, then P (s 0 = 1) = 1 and P (s 0 = 2) = 0.
We now re-write Eq. 21 as  P (s 3 |s 2 )P (n 2 |s 2 )× s1 P (s 2 |s 1 )P (n 1 |s 1 ) s0 P (s 1 |s 0 )P (s 0 ) ... . (22) This equation represents the probability of observing n i products in a given time level i, P (n i |s i ) which depends on the state of the catalytic site at that time, and the probability of switching from s i−1 to s i , P (s i |s i−1 ). On account of the Markov assumption, the probability of observing a state at each time level only depends on its realization at the previous time level.
Hence, for a given initial condition at time t = 0 (i.e. P (s 0 )) and the trajectory of products formed generated at all time levels (n 1:M ), we can calculate the probability of observing a trajectory n 1:M by summing over s i 's starting from i = 0 and moving up to i = 1, then i = 2, 3, .., M − 1 in a recursive fashion and then, finally, over s M . This protocol is called Forward Filter [7], and is summarized below P (n 1 , s 1 ) = P (n 1 |s 1 ) s0 P (s 1 |s 0 )P (s 0 ) P (n 1:2 , s 2 ) = P (n 2 |s 2 ) s1 P (s 2 |s 1 ) × P (n 1 , s 1 ) In the section above, we devised a method for calculating the probability of observing a given trajectory comprised of a number of product formation events observed at M time levels, Eq. 22. Using this expression, we can calculate the information contained in those trajectories where P (n 1:M ) is the probability of the trajectories (that we sample). Here, the summation is performed over all the possible trajectories which can be observed in the model (see main text). Summing over this large space is computationally expensive. This can be avoided by calculating the information over a subset of L sampled trajectories, C(n 1:M ), as shown Since trajectories in C n 1:M are observed in the stochastic simulation describing the same model, these will be sampled from the same distribution, P (n 1:M ) as our model. For each of the trajectory sampled stochastically, the probability of observation of that trajectory can be calculated using the method described in the section above Eq. 22. Finally, the information over these subset of trajectories, I is then calculated using Eq. 25. Therefore, for large values of L, we expect, within acceptable error, that lim L→inf I = I.
Since only few trajectories have appreciable probability, the condition for convergence as shown by Eq. 26 can be satisfied for a reasonable value of L. For parameters λ 1 = 60, λ 2 = 100, P 1→2 = P 2→1 = 0.04, we observed that the information converges around L ≈ 150, see Fig. 4.
Moreover, we also calculated the ratio of information contained in the model with the length of the trajectories sampled (i.e., the number of time levels observed in the trajectory). The calculations were performed over 800 trajectories sampled from stochastic simulation (such that the information calculated is converged) and we observed that the ratio does not change significantly with the length as shown in Fig. 5.
These results ensures that the stochastic simulation is sampling from a stationary Markov process (on account of the steady state assumption) and that 150 trajectories are sufficient for the information to converge.