The model
We show in Fig. 1a scheme of the model considered in the paper. A constant external stimulus can either elicit a single pulse of TF nuclear traslocation (amplitude encoding) or a sequence of pulses (frequency encoding). Namely, during the simulation time, \(0\le t\le T=100\, {\text{min}}\), the TF nuclear concentration, [TF], is given by:
(1)
$$\begin{aligned} [TF](t)&= \zeta (t)+\sum _i 100\, H(t-t_i)H(t_i+1\,{\text{min}}-t), \, \, \, \, \mathrm{in \, frequency \, encoding,} \end{aligned}$$
(2)
where H is the step function (\(H(x)=1\) for \(x\ge 0\) and 0 otherwise) and \(\zeta\) and the inter-pulse time interval, \(\tau \equiv t_{i+1}-t_i\), are random variables described later in more detail. The external stimulus strength, \(I_{ext}\), determines \(A_{TF}\) and \(\langle \tau \rangle\) for each encoding, respectively, according to:
$$\begin{aligned} A_{TF}&=100\frac{I_{ext}^{h}}{I_{ext}^{h}+EC_{50}^{h}} , \end{aligned}$$
(3)
$$\begin{aligned} \langle \tau \rangle&=T_{{\text{min}}}\left( 1+\kappa \exp \left( -bI_{ext}\right) \right) , \end{aligned}$$
(4)
where \(EC_{50}\), h, \(T_{{\text{min}}}\), \(\kappa\) and b are fixed parameters. In the model all the concentrations, \(I_{ext}\), and the parameters in Eqs. (3) and (4), with the exception of \(T_{{\text{min}}}\), are dimensionless. The rationale for choosing Eqs. (1)–(4) as well as the meaning of the various parameters are explained later in this Section. We have not included a detailed dynamical system to derive TF from \(I_{ext}\) since it would only add to the transient of the TF dynamics without affecting the (100 min) mRNA time integral that we take as the output of the process ( see e.g., the example of Msn2 in yeast where the transients last for \(\sim 5-8\,{\text{min}}\)7).
[TF] modulates the transcription rate according to the model depicted in the second row of Fig. 113,16, where \(P_0\) and \(P_1\) represent the two (effective) states of the promoter and transcription proceeds only when the promoter is in the \(P_1\) state at rate, \(k_2 [TF]^nP_1/(K_d^n+ [TF]^n)\) (as explained later, we also use \(P_0\) and \(P_1\) to denote the probabilities that the system is in one or the other state), while mRNA is degraded at rate, \(d_2\). This part of the model as well as the indicator of gene expression (the accumulated mRNA) are the same as those used in our previous study16 where we performed an extensive search of optimal parameter values for information transmission through the transcription step. The parameters of the transcription step used in the present paper were chosen based on this previous study. We describe the model in more detail in what follows.
The model. An external stimulus elicits a single pulse or a sequence of pulses of nuclear TF translocation in amplitude and frequency encoding, respectively. The stimulus strength, \(I_{ext}\), (the input) determines the mean concentration, \(A_{TF}\), of the single TF pulse and the mean time, \(\langle \tau \rangle\), between successive pulses in amplitude and frequency encoding, respectively, according to the expressions inside the blue box (Eqs. 3 and 4). The rationale for these choices as well as the meaning of the various parameters are explained later in Methods. The resulting nuclear TF concentration, [TF], depicted for each encoding in the blue box, modulates the transcription rate according to the two-state promoter model presented in the second row which, in turn, determines the mRNA time integral that is taken as the indicator of gene expression, i.e., the output (third row). The transcription model is the same as the one considered in our previous work16. We recall here that \(P_0\) and \(P_1\) denote the two promoter states as well as the probability that the promoter is in each of these states, that the two arrows connecting \(P_0\) and \(P_1\) in the scheme represent the transition between the two states while the arrow that goes from \(P_1\) to mRNA is used to mean that the mRNA production proceeds at a rate given by \(k_2 [TF]^nP_1/(K_d^n+ [TF]^n)\).
Transcription
Transcription is modeled equally for frequency and amplitude encoding using the effective two-state promoter model13,16 depicted in the second row of Fig. 1. Using \(P_0\) and \(P_1\) to denote the probabilities that the promoter is in one or the other of its two states, noticing that \(P_0=1-P_1\), and recalling that the arrow connecting \(P_1\) with mRNA in the scheme means that the mRNA production proceeds at rate, \(k_2 [TF]^nP_1/(K_d^n+ [TF]^n)\), the dynamical equations that are used to simulate the model read:
$$\begin{aligned}&\dot{P}_1=\frac{k_1 [TF]^n}{K_d^n+[TF]^n}-\left( \frac{k_1 [TF]^n}{K_d^n+[TF]^n}+d_1\right) P_1,\end{aligned}$$
(5)
$$\begin{aligned}&X=N\xrightarrow []{\frac{{k_2} [TF]^nP_1}{K_d^n+ [TF]^n}}{X=N+1},\, {X=N}\xrightarrow []{d_2 X}X=N-1 , \end{aligned}$$
(6)
where [TF] is the nuclear TF concentration (in dimensionless units), \(P_1\) represents the probability that the promoter is in its active TF-bound state and X(t) is the random number of mRNA molecules at time, t, that can take on the values \(N=\{0,1,2\ldots \}\). Equation (5) is a deterministic equation for the probabilities and Eq. (6) is a stochastic equation for the random variable, X. The extensive search of “optimal” parameter values performed before for the transcription model16 showed that MI achieves its maximum in the same bulk region of parameter space for amplitude and frequency encoding. The search did not yield sharp maxima either16. Thus, in this paper we use a fixed set of parameter values for the transcription step within this previously determined region, limiting the exploration of parameters to those determining the mapping from the external stimulus or the \(I_{ext}\) distribution. The simulations are then performed with time step \(dt=0.1\,{\text{min}}\), integration time \(T=100\,{\text{min}}\), \(n=10\), \(K_d=40\) (in the same dimensionless units as [TF]), \(k_1=1/{\text{min}}\), \(d_1=0.01/{\text{min}}\), \(k_2=10/{\text{min}}\) and \(d_2=0.12/{\text{min}}\) and using Euler’s method to integrate Eq. (5). In all cases, the output is:
$$\begin{aligned} O= \int _0^T dt\, X(t) . \end{aligned}$$
(7)
From the external stimulus to the nuclear TF
For amplitude modulation, [TF] is modeled as a single pulse of 10 min duration and (dimensionless) amplitude (see Eq. (1)):
$$\begin{aligned} A=A_{TF}(I_{ext})+\zeta , \end{aligned}$$
(8)
with \(\zeta\) a Gaussian distributed random variable of standard deviation, \(\sigma _\zeta =10\), and \(A_{TF}\) related to the external input strength, \(I_{ext}\), (typically, the dimensionless concentration of an external ligand), by the cooperative Hill function (3). Equation (3) commonly describes the input-output patterns observed in signaling cascades46 and, in the context of the present model, aggregates in one step the various processes that go from the external stimulus to the TF’s nuclear fraction. \(EC_{50}\) and \(I_{ext}\) are measured in the same dimensionless units, which not necessarily coincide with those of [TF].
For frequency modulation, [TF] is modeled as a sequence of 1 min-duration square pulses of (dimensionless) amplitude 100 plus a random variable, \(\zeta\), as in the case of amplitude encoding (see Eq. (2)), and stochastic interpulse time intervals,
$$\begin{aligned} \tau = T_{{\text{min}}}+\eta , \end{aligned}$$
(9)
with \(T_{{\text{min}}}\) fixed31 and \(\eta\) exponentially distributed with rate parameter exponentially dependent on \(I_{ext}\)42,45 so that:
$$\begin{aligned} T_{IP}(I_{ext})&\equiv \langle \tau \rangle =T_{{\text{min}}}\left( 1+\kappa \exp \left( -bI_{ext}\right) \right) , \end{aligned}$$
(10)
with \(\kappa , b >0\). In this case b is measured in the inverse of the (dimensionless) units of \(I_{ext}\).
Input, output and mutual information
In this paper we compute MI between the external stimulus strength, \(I_{ext}\), and the mRNA produced, O, over the time course of the simulation (Eq. (7) with \(T=100\,{\text{min}}\)). MI can be written in terms of the marginal and joint probability densities of these two variables, \(p_I\), \(p_O\) and \(p_{I,O}\), respectively, as:
$$\begin{aligned} MI=\int p_{I,O}(I_{ext},O)\log \left( \frac{p_{I,O}(I_{ext},O)}{p_I(I_{ext})p_O(O)}\right) . \end{aligned}$$
(11)
MI is computed numerically using an \(I_{ext}\) distribution of compact support, \([I_m,I_M]\), and the Jackknife method47,48. This method is commonly used to derive entropy estimates since it reduces the bias of the Maximum Likelihood (ML) estimator of \(x\log {x}\)47,48. It consists of taking a linear regression of the ML estimator for various sample sizes, with the slope being the bias of the ML estimator and the intercept being the Jackknife estimator. For the computation we used \(T_{{\text{min}}}=5, 10\, {\text{min}}\), various values of \(\kappa\) and the \(I_{ext}\) distribution, \(p_I\), derived from either one of these expressions:
$$\begin{aligned} I_{ext}&= x,\end{aligned}$$
(12)
$$\begin{aligned} I_{ext}&= \exp \left( 4(x-0.5)\right) , \end{aligned}$$
(13)
with \(0\le x\le 1\) a stochastic variable distributed according to the Beta-distribution:
$$\begin{aligned} f(x)&= \frac{x^{\alpha -1}(1-x)^{\beta -1}}{\textrm{B} (\alpha ,\beta )},\quad 0\le x\le 1, \end{aligned}$$
(14)
for different choices of \(\alpha , \beta >1\) s.t. \(\alpha + \beta \ge 3\), so as to obtain different values for the median of \(I_{ext}\). The Beta-distribution (Eq. 14) is often used to model fractional quantities. It is the posterior distribution that is obtained when applying a Bayesian approach to estimate the probability of a Bernoulli process from \(\alpha +\beta -2\) observations49. Assuming that the external stimulus corresponds to a ligand that, upon binding to a membrane receptor, triggers a signaling cascade, the probability of the Bernoulli process can be interpreted as that of a ligand’s molecule being within a certain (interaction) distance from the membrane (when Eq. (12) is used). This probability will be proportional to the ligand’s concentration. Fixing \(\alpha +\beta\) and choosing different values of \(\alpha\) will then correspond to having different (mean) ligand concentrations nearby the cell. This is illustrated in Fig. 2.
Cumulative distribution function (CDF) of the stimulus strength, \(I_{ext}\), obtained with Eq. (12) (a), (c) and Eq. (13) (b) combined with Eq. (14) for 10 equally spaced \(\alpha\) values increasing from left to right and \(\alpha + \beta\) fixed ((a), (b): \(1.2\le \alpha \le 3\) and \(\alpha +\beta =4\); (c): \(1.4\le \alpha \le 5\) and \(\alpha +\beta =6\)). Varying \(\alpha\) and \(\beta\) for increasing values of \(\alpha +\beta \ge 3\) yields broader ranges of the \(I_{ext}\) median and standard deviation (with smaller values for the latter: \(\sim 0.19-0.22\) in (a) and \(\sim 0.14-0.19\) in (b)). Values such that \(\alpha +\beta = 4,5\) present a good balance of medians that cover relatively well the [0,1] interval with relatively invariant and not too small variances.
As shown in Fig. 2, the range of median \(I_{ext}\) values that can be explored by changing \(\alpha\) when Eq. (12) is used varies approximately linearly between 0 and 1. In order to explore a wider range of values, we decided to probe Eq. (13) for which the range varies logarithmically (Fig. 2b). Increasing values of \(\alpha +\beta \ge 3\) yield more widely spread medians and smaller values of the standard deviation as illustrated in Fig. 2. In the paper we show the results obtained with the distributions of Fig. 2a because their medians cover reasonably well the whole interval, [0, 1], with a relatively invariant and not too small standard deviation (between 0.19 and 0.22).
Computation of the distribution function of the intermediaries of the response
In the approach followed in this paper the intermediaries of the response are the amplitude, A, or the interpulse time, \(\tau\), (or, equivalently, the frequency \(\nu \equiv 1/\tau\)) of the TF’s nuclear concentration in the case of amplitude and frequency encoding, respectively. The amplitude, A, is given by Eqs. (3), (8) where \(\zeta\) is normally distributed with standard deviation, \(\sigma _\zeta\). This implies that A is Gaussian distributed with mean \(A_{TF}\) and standard deviation, \(\sigma _\zeta\). The interpulse time, \(\tau\), is given by Eq. (9) with \(\eta\) exponentially distributed with mean, \(\langle \eta \rangle = T_{IP}-T_{{\text{min}}}\), where \(T_{IP}\) is given by Eq. (10). A and \(\tau\) (or \(\nu\)) can be thought of as the “input” of the transcription model. Given that their distributions are known for a given value of the external stimulus, \(I_{ext}\), and that the \(I_{ext}\) distribution is known as well, it is possible to derive the A and \(\tau\) or \(\nu\) distributions that are fed into the transcription model in our approach. In particular, in the paper we show the results obtained for the choice of \(I_{ext}\) distribution given by Eqs. (12) and (14). In such a case, the cumulative distribution functions (CDF) of A and \(\nu\) can be computed as:
$$\begin{aligned} CDF_{A}(A) = \int ^A_0 da\int _0^1 dI_{ext} f(I_{ext}) \frac{e^{-\frac{\left( a-A_{TF}(I_{ext})\right) ^2}{2\sigma _\zeta ^2}}}{\sqrt{2\pi \sigma _\zeta }}, \end{aligned}$$
(15)
$$\begin{aligned} CDF_{\nu }(\nu ) = 1-\int ^{1/\nu }_{T_{{\text{min}}}} dt \int _0^1 dI_{ext} f(I_{ext}) \frac{ e^{-\frac{t-T_{{\text{min}}}}{T_{IP}(I_{ext})-T_{{\text{min}}}}}}{T_{IP}(I_{ext})-T_{{\text{min}}}}, \end{aligned}$$
(16)
with f given by Eq. (14), \(A_{TF}\) by Eq. (3), \(T_{IP}\) by Eq. (10) and where we are assuming in Eq. (15) that the integral in a from \(-\infty\) to 0 is negligible for most values of A. We use these expressions to compute numerically these two CDFs (which are shown in Fig. 4c,d in the “Results” section).