Maximization of information transmission influences selection of native phosphorelay architectures

Mechanistic models

PR shown in Fig. 1 represent more than 99% of all PR circuits present in the fully sequenced genomes of more than 9,000 organisms. This information was obtained from Salvado et al. (2015).

We used a mass action description of reactions to create a mechanistically more detailed model for each architecture (Alves et al., 2008). These models are given in Appendix S1. They were parameterized with experimentally determined parameter values and protein concentrations compiled from the primary literature. We used these models to create an atlas of performance with information transmission through the PR as a function of protein amount, parameter values, and circuit architecture. The Mathematica notebooks implementing the models and analysis are provided in Supplemental File S5.

Estimating protein concentrations

Proteins occupy approximately 20% of cell volume across the tree of life (Dill, Ghosh & Schmit, 2011; Milo & Phillips, 2015). Taking into account average protein sizes, average protein masses and factoring in cell volumes one can estimate that total protein concentration in cells is of the order of ∼1 µM (ranging between 0.4−1.4) (Uri & Milo, 2010), or 6. 23 × 10⁶ proteins per µm³. We then use the average cell volume for the different cell types to estimate the total amount of proteins per cell.

In order to estimate biologically relevant protein amounts for each circuit topology we used whole-proteome relative abundance determinations (Table S1) reported by PaxDB (Wang et al., 2015). The database information is given in parts per million. By multiplying these number by total number of cell proteins, we calculate how many proteins exist for each experimental PR system, within the experimental error.

With this information we can further estimate typical orders of magnitude for the ratios between protein abundances within a PR circuit. In general, we find that protein abundances and ratios between protein abundances have a limited range in PR (Table S1).

Parameter values for the mass action models

To find experimentally determined parameter values for the individual reactions of the PR, we searched Medline and the Biomodels database (Glont et al., 2018). We searched for previous mathematical models for TCS and signal transduction PR. In addition we also search the primary literature for known circuits that have been characterized biochemically and collected all different parameter values for the various reactions in the PR circuit. This information is compiled in Table S2, revealing that kinetic parameters for corresponding reactions in experimentally well-characterized systems are quite similar and, in most cases, within the same order of magnitude.

Steady state concentrations and stability

We obtained steady state concentrations by numerically solving the models for each set of parameter values. We then calculated the Jacobian of the ODE system and the eigenvalues of that Jacobian to determine the stability of the steady state. We compared the steady state stability of the PR circuits in two ways. As we selected for sets of parameter values that generate systems with stable steady states (see methods), all real parts of the eigenvalues of the steady state are negative. Taking this into account we compared the minimum of the real parts of the eigenvalues in each pair of models. This comparison allows us to compare the fastest time scale in which the two circuits respond to a transient perturbation to their steady states. We also compared the maximum of the real parts of the eigenvalues between the two circuits in each pair. This allows us to compare the slowest time scale in which the two circuits respond to a transient perturbation to their steady states.

Logarithmic gains and parameter sensitivities

To estimate the response of a physiological variable to changes in environmental signals or model parameters, we compute (Alves et al., 2008) ${\displaystyle L\left(X,\theta \right)=\frac{dLog\left(X\right)}{dLog\left(\theta \right)}={\left(\frac{\theta }{X}\times \frac{dX}{d\theta }\right)}_{SS}}$

Here, ss indicates evaluation at a reference steady state.

Signal amplification

To estimate signal amplification we calculated the logarithmic gains of the concentration of the phosphorylated form of the final response regulator (RR2) with respect to the cognate signal of the system (Signal), as given by: ${\displaystyle L\left(RR2P,Signal\right)=\frac{dLog\left(RR2P\right)}{dLog\left(Signal\right)}=\frac{Signal}{RR2P}\frac{dRR2P}{dSignal}}$

In our case the Signals are the rate parameters for the self-phosphorylation and self-dephosphorylation reactions of the SK.

Sensitivity to the total amount of protein

To estimate sensitivity of the different phosphorylated forms of the PR proteins (Prot-P) with respect to the total amount of circuit proteins (Pr_tot), we calculated the logarithmic sensitivities, as given by: ${\displaystyle L\left(Prot-P,P{r}_{tot}\right)=\frac{dLog\left(Prot-P\right)}{dLog\left(P{r}_{tot}\right)}=\frac{P{r}_{tot}}{Prot-P}\frac{dProt-P}{dP{r}_{tot}}}$

Global sensitivity to fluctuations in parameters

To estimate sensitivity of the different phosphorylated forms of the PR proteins (Prot-P) with respect to fluctuations in each of the parameters (Par _i), we calculated the logarithmic sensitivities, as given by: ${\displaystyle S\left(Prot-P,Pa{r}_i\right)=\frac{dLog\left(Prot-P\right)}{dLog\left(Pa{r}_i\right)}=\frac{Pa{r}_i}{Prot-P}\frac{dProt-P}{dPa{r}_i}}$

Then, for each protein, we calculate the norm of the vector whose coordinates are the individual sensitivities of that protein, as described in (Savageau, 1971a; Kitano, 2002). Thus, for architecture Mi, the aggregated sensitivity is given by ${\displaystyle S_{aggregate,Mi}=\sqrt{\sum _{i=1}^{n}S{\left(Prot-P,Pa{r}_i\right)}^2}}$

This is an aggregated measure of the robustness of circuits to fluctuations in parameter values (Alves et al., 2008).

Response times

We performed two independent experiments per architecture and per set of parameter values to calculate the response times. First, we start with the system fully dephosphorylated and run a time course simulation using the same scan procedure for the self-phosphorylation and self-dephosphorylation rate constant of SK. We then calculate the time it takes each of the architectures to reach 90% of the new steady state values. We repeated this experiment starting with fully phosphorylated proteins.

Amount of information transmitted through the circuit

To calculate the information transmitted through the PR we simulated that system to steady state using Gillespie’s algorithm for stochastic simulation one hundred times per set of parameter values (Gillespie, Hellander & Petzold, 2013). The proxy for changes in the environmental information was taken to be the number of phosphorylated molecules of the SK domain, as this is the sensing domain of the PR. The output of the PR was considered to be the number of phosphorylated final response regulators. A high correlation between the variations in the amount of phosphorylated SK and RR2 implies that information is transmitted in higher amounts through the circuit.

To calculate the amount of information being transmitted through the circuit we started by running twenty five stochastic simulations with each set of parameter values and initial conditions. At steady state, we sampled each simulation at one second intervals, running the steady state for 20 min. We pooled together the simulations and calculated the individual relative frequencies, for the number of SKP, RR1P, HptP and RR2P molecules. We also calculated the paired relative frequencies for (SKP, RR2P). We then used these numbers to estimate the information being transmitted through the circuit by applying the standard formula for mutual information between the two variables (Bekker, Teixeira de Mattos & Hellingwerf, 2006; Lan & Tu, 2016): ${\displaystyle MI\left(SKP\to RR2P\right)=\sum _{SK{P}_i}\sum _{RR2P_j}p\left(SKPi,RR2Pj\right)log\left(\frac{p\left(SKPi,RR2Pj\right)}{p\left(SKPi\right)p\left(RR2Pj\right)}\right)}$

In these formulas, p(x) is the relative frequency of x, and p(x, y) is the joint relative frequency of x and y. $MI\left(SKP\to RR2P\right)$ estimates the amount of information being transmitted through the circuit.

We calculated $MI\left(SKP\to RR2P\right)$ as a simultaneous function of the rate constants for the phosphorylation and dephosphorylation reactions of the SK and RR2. To present this information in a way that is visually easier to interpret, we calculate an aggregated mutual information defined as: ${\displaystyle M{I}_{agregated}\left(SKP\to RR2P\right)=\sum _{k=kmin}^{kmax}MI\left(SKP\to RR2P\right)}$

Here, k represents the rate constant that is not explicitly represented in the plots.

Numerical comparison of mechanistic models

We took a reference architecture and generated a latin hypercube sampling for the parameter values of the mechanistic models of architecture M1. We approximated the experimentally determined protein abundances to the order of magnitude and then also generated a latin hypercube sampling strategy, which we combined with that for the parameter values. Overall, we performed several hundred million sets of simulations.

Mathematical modelling of the Spo0 and Sln1 phosphorelays

Detailed mechanistic models were created for the Spo0 and Sln1 phosphorelays, as described in Appendix S1. These models use experimental data from Narula et al. (2015) for Spo0 and from the SGD (Hirschman, 2006) for Sln1.

Then, for each of the PR, four additional equivalent models were created, each assuming that the PR would have a different architecture. The parameter values for these alternative architectures were considered to be the same as those for the native circuit architecture. The concentrations for the proteins were optimized by first simulating the original model at different signal intensities and calculating the steady state phosphorylation state of the final RR. Using these results we then performed the same experiments for each alternative architecture and allowed for the proteins that got fused or separated with respect to the cognate architecture to change in concentration by up to one order of magnitude below or above that in the natural architecture. We then selected the concentrations that led to the most similar signal-response curves using a least square minimum criteria for the differences between phosphorylated final RR.

Calculating the cost of alternative architectures for the Spo0 and Sln1 phosphorelays

To calculate the cost of alternative architectures for the Spo0 and Sln1 we first counted the number of amino acids in the sequence for each protein in the original architecture. Then, we assumed that any alternative architecture would result from fusing or splitting the original proteins while conserving the same number of amino acids. Finally, we multiplied the number of amino acids of the individual proteins in the phosphorelay by the number of proteins in the cell and used this number as a proxy for the cost of each architecture. This makes it so that the cost of maintaining one copy of the PR circuit is the same between the alternative architectures. Thus, differences in the metabolic cost for alternative architectures of the Spo0 (or Sln1) PR circuit are simplified to the differences in ATP consumption rate. As we also adjusted concentrations of each alternative architecture to make the signal response curves of the alternative architectures as similar as possible to that of the native architecture (see below), this makes differences in the cost of circuit maintenance equate to differences in protein concentration.

Software

All models and analysis were done using Mathematica (Wolfram Research I, 2014). The notebooks containing all code to generate each figure are given as supplementary data pack S1.

A performance atlas for information transmission in the parameter space of PR architectures

With the hypothesis that information transmission is important in the selection of PR architecture, we created a data-driven atlas that describes how architecture, protein abundance, and parameter values influence information transmission in PRs. To create that atlas we compiled experimentally determined parameter values and proteins abundances for as many PR systems as we could find in the primary literature and in public databases (see methods and Appendix S1). Then we used those parameters and protein abundances to calculate the capacity to transmit information of each alternative architecture with all possible combinations of protein abundances and kinetic values. Finally we ranked the five architectures with respect to increasing capacity of information transmission. The atlas is summarized in Table S3.

Experimentally determined protein abundances are consistent with amount of transmitted information being an important driver of PR architecture selection

We identified and obtained protein abundances for seventeen PR systems that are present in the PAXDB database (Wang et al., 2015) and belong to seven different organisms (Table S1). Then we matched these abundances with the performance in the atlas of amount of transmitted information as a function of parameter values, protein abundance and circuit architecture given in Table S3. In 16 out of the 17 cases we find that the native architecture of the PR systems is consistent with maximizing the through-circuit transmitted information under experimentally known parameter conditions (Table 3).