Self-organization of signal transduction

Example system: GPCR signaling in a submembrane compartment

Figure 1A shows the example system, a submembrane compartment with a GPCR (G-protein coupled receptor) pathway from a mouse embryonic fibroblast, with ISO as input (extracellular ligand to β(2)-adrenergic receptors) and the phosphorylation of a protein VASP as output. This model was implemented as an ODE model with 23 reactions and 27 molecular species, derived from 12 initial concentrations (cf. Table 1, Table 2). The parameters were adapted to experimental biological data (not shown¹¹). In this subsystem, the central cAMP response often follows a plateau curve, i.e. a rise to steady-state, but cAMP transients which are typically observed in cytoplasm may also occur¹¹. The dose-response transfer functions were derived as in⁶.

Objective function

A biological signal transduction system is defined by its state variable vector x, the set of all kinetic rate and initial concentration parameters.

We hypothesize that an efficient signal transmission would maximize the response coefficient R_C,S (the response of species C to input S) defined as

with concentration change of target C and input S from baseline (t=0) to signal time t. For R_C,S, values < 1 show signal loss, with 1 for perfect transmission, and values > 1 showing signal amplification. We may also optimize for the slope s of the sigmoid at half-maximum concentration. This is equivalent to maximizing R_C,S, provided that the input signal remains entirely between the upper and lower boundaries of the sigmoid (Figure 1C). By optimizing for s, additional to R_C,S, we may force the system to implement a switch-like function instead of a more linear function. However, shifting the sigmoid function to the left or to the right is more important as the slope in our models.

In addition, we optimize for reaction time (delay to steady-state). Steady-state is defined, pragmatically, as relative change of less than 2% over 100s. The delay (d_S) is computed for 90% (EC90) of steady-state. We may now define an aggregate objective function:

to select the system state variables that minimize delay and maximize response.

In addition to signal transmission from extracellular concentration changes onto steady-state concentrations, such as they typically occur for temporally integrating proteins like transcription factors, we also look at signal transmission by transients, i.e. peak concentration, in response to extracellular signals. In this case, we minimize the delay to peak value, maximize the response at peak value, and minimize the response at steady-state value.

Delay vs. efficiency trade-off

The computation of signal transmission efficiency will be explained here for a single input-output pathway of cAMP/PKA-mediated transmission in a cellular membrane compartment. It is clear that a complex signaling system may have several inputs, and a large number of outputs or target proteins, and this is especially the case for cAMP-mediated signaling. Nonetheless we will focus on the simple case here to explain the basic principle. The input is a membrane receptor, a G-protein coupled receptor (GPCR), β(2)AR, which is activated by an extracellular pharmacological agonist ISO (isoproterenol); the output is a membrane protein, VASP, which promotes actin filament elongation. Activation of the receptor selectively inhibits VASP by phosphorylating VASP to pVASP via protein kinase A (PKA) activation. In Figure 1B the original ISO/pVASP transfer function is shown, which we use here for further optimization. The system was trained for a signal distribution between 10nM and 1 µM ISO adjusting both kinetic rate and concentration parameters (Figure 1A). The optimization method used is a simplex algorithm¹². Figure 1B shows the transfer function before and after adjustment, maximizing for R_C,S, d_S or the combined function f. The results are summarized in Table 3. We see that optimizing for the response coefficient alone shifts the function to the right to better cover the input range. Optimizing for the delay alone shifts it to the left, speeding up signals in the lower range (which are slower). Both objectives together (equally weighted) are fairly close to the original, biologically validated curve, with an improved R_C,S.

In general, biochemical reactions are faster at higher substrate concentrations, but the relative concentration change in response to an increase in the enzyme or the binding partner is less. This is a fundamental trade-off between delay and transmission efficiency that may define an optimal operating range for a signaling system and be of relevance in disease processes⁶. The results obtained with this experiment are simple, intuitive and encourage continuing to explore the basic idea.

Evolutionary learning of kinetic rate parameters

In principle, we may use all parameters in a system, concentrations or kinetic rates, to maximize signal transmission. But the evolution of protein structure and interactions shows that it is fine-tuning of molecular kinetic parameters which is subject to evolutionary learning, while concentrations are often regulated adaptively in each cell. Mass-action kinetics approximate molecular kinetic parameters, even though there are significant sources of uncertainty, such as the stochastics of molecular interactions. In the following, we explore the idea that kinetic rate learning operates on evolutionary time-scales, and that biologically attested signal transduction pathways contain reaction rates which are optimal in terms of signal transmission efficiency. We use known concentration ranges, specific by cell type, together with kinetic rate optimization.

To explore the parameter space, we drew 1662 values for all kinetic rate parameters (kon, koff, kcat) from a distribution of 20% to 500% of the original parameters, and calculated R_C,S and d_S for the corresponding models, relative to an improved signal distribution from 1nM to 1µM (Figure 2A). We find that there are parameter combinations which greatly improve efficiency of the signal transmission function. The basic distribution of a uniform low value of R_C,S for fast d_S and a wide variability in R_C,S above a certain threshold in d_S is robust against different types of signaling input (cf. also Figure 3A). This may, however, be highly dependent on the reaction network that underlies the transfer function. We have not further explored this question. To test for robustness of these systems against variability of concentration, we repeated experiments for 100 systems with 20% variation of original concentrations, which corresponds to generally accepted noise levels (cf.^13,14). As expected, this low variation did not significantly affect the quality of a set of kinetic rates (Supplementary table).

Figure 2.

A. Distribution of systems according to R and d values. ~ 1600 different systems states were randomly generated with k-parameter variations (20–500%), R_C,S and d values were measured with the signal set shown on top. B. Kd and KM parameters for high efficiency (red, f < 2) and low efficiency (green, f > 5) systems. C. Cross-correlation of Kd values.

Figure 3.

A. (top) Dynamics of selected systems with low (f_transient < 3, blue) or high (f_transient > 10, red) propensity for a transient response. (bottom) Distribution of concentration parameters according to d^p, $R_{\mathit{\text{C,S}}}^p$ (grey) or R_C,S, d (red). B. Concentration changes in response to optimizing for transients, with plateau signaling (left) and strong transients (right), sorted by f_transient. C. Concentration shifts (orange=high, green=low) in the biochemical reaction network for a high transient system, as averaged from B. Conforming to intuition, we see that the earlier, driver complex (G-protein) is high for transients, which increase the on-slope for cAMP production.

We further analysed the parameter combinations with different signal transmission efficiency. In Figure 2B and C, we distinguish low and high efficiency signal transmission. Interestingly, we find, with respect to signal efficiency of the transfer function, that all parameters are “sloppy” (allow a wide variation), there are no “stiff” (low variation) parameters¹⁵. Standard deviation for all parameters in our case is similar (Figure 2B). Since it has been argued that optimization to experimental data yields reactions which allow more variance than others, as an indication of their influence on the signal transmission pathway, this analysis seems to contradict this effect. Possibly, these results pertain mostly to parameter variation that results from matching a networked system with many species to selected time-series data for only few species, which may behave differently from general optimization.

Co-regulation of protein concentration as an adaptive response to perturbation

Co-regulation of protein expression in cellular systems is important in disease progression and often a problem in targeted interventions. Here we are exploring the question of self-organization of protein concentration after a perturbation that reduces one protein to only 10% of its previous concentration. Keeping kinetic rates fixed, all concentrations in the system are allowed to adjust until optimality of signal transmission and delay is reinstated. There is a number of interesting observations here (cf. Figure 4A), which relate to the biological reality. For instance, reducing PDE4B causes much regulation in other proteins, but it is almost never targeted. In contrast, reducing PP1 has little effect, but PP1 is frequently responsive to other proteins. Reducing PKA, RGS and AC6 leads to widespread down-regulations, to maintain sensitivity of signaling, but reducing the receptor beta-2 leads to up-regulation.

Figure 4.

A. Concentration changes in response to total protein concentration reduction. Selected concentrations were reduced to 10% (shown on the x-axis). Learning was applied until f values were improved (shown on top for each experiment). Colors show the relative adjustment of all concentrations on the y-axis. B. Transfer function response to shifts in input signal range (shown on top). C. Concentrations changes in response to extracellular signal shift. Matrix is constructed as in A.

There are many individual adaptations which can be interpreted to maintain the sensitivity and responsivity of the small molecule cAMP. The results show that protein regulation is highly sensitive to positions and roles of individual proteins in the reaction network in transmitting the signal. This problem may also be amenable to a more principled mathematical analysis¹⁶. Another idea would be to rank reaction systems defined by kinetic parameters as in Figure 2 by how well they adapt to perturbations.

In our example, the quality of the readjustment is not always the same, but the simple optimization scheme that we use may easily produce suboptimal solutions (local minima). Concentration changes may be caused by genetic up- or down-regulation, secretion and re-uptake, increased degradation, RNA interference, etc. and are therefore not easy to model from the standpoint of mechanistic biological modeling, which would need modules for all biologically attested processes. A unified perspective by a set of constraints and a set of objectives, such as has been envisaged here, may lead to better predictive results and may also be used as a guiding principle in constructing mechanistic models. Since there are intricate biological processes of adjusting concentrations, we may assume that in the cell optimal solutions are more easily found.

Optimal signal transmission depends on the signaling level

From the standpoint of disease modeling, an unusual protein concentration may be an adaptive response where a still functioning cell in a dysfunctional external signaling environment struggles to keep signal transmission efficient to support cellular function. In such a case, targeting this protein by pharmacological intervention will lead to co-regulation on other proteins. In general, as well as in real biological signaling systems, it is the localization of the input signaling range and the distribution of signals that the system transmits which are important for optimization. If we select the signal set that we optimize for in the right way, the input range will move towards the loglinear range, i.e. between the lower and upper input boundaries (Figure 4B). As a result, a number of internal concentrations will change. Figure 4C shows the relative concentration changes that result from a shift in input range. Interestingly, the low shift requires a strong reduction in RGS, and upregulation for PDE4B, among other effects. GsaGDP (the activating G protein) is strongly increased, and GiGDP (the inhibitory G protein) is decreased. In the high shift, GiGDP and PP1 (the phosphatase which decreases pVASP) are stronger, and the beta-2 receptor density is much increased. In the future, it will be interesting not only to calculate these effects based on different optimization measures, and in multiple input-output scenarios, but also to compare this with attested cases of biological adaptation. By reverse engineering, we may infer an optimization measure from a sufficient number of attested co-regulations. For instance in addiction, protein co-regulation as a form of sensitization is well-attested¹⁷, but also in cancer where intercellular signaling (e.g. by cytokines¹⁸) is affected. Gene expression data may then be mined not only for evidence of the mechanics of genetic regulatory pathways but also for evidence of shifts in extracellular signaling, which cause altered protein expression.

Role of transients

The appearance of a transient vs. a plateau signal (or even a dampened oscillation) in response to sustained signaling depends on the construction of the biochemical reaction network (negative feedback interactions) and its parameters. We show that we can also train a system for the appearance of a transient response to a sustained signal. This means to search for high response at peak, low delay to peak, but also a low response at steady-state (see Section “Objective function” and Figure 3A). This requirement of invariance for the steady-state is sometimes considered a form of “robustness” or “homeostatic regulation” of the cellular response^19,20, and signaling by transients is widely regarded as an important mechanism. Here we found that concentration adjustment is quite sufficient to acquire a switch from transient to plateau response, with a high variability of response shape dependent on concentration, but also on the size of the signal (Figure 3A). In Figure 3B we show the distribution of concentrations that achieve high or low propensity for transients, as indicated by f_transient. We focus on concentrations with fairly uniform up- or down-regulation to create an experimental graph in Figure 3C. The simplicity of Figure 3C undermines the notion, as proven by the random search optimization, that deviations from up- or down-regulation for individual concentration may not be irrelevant noise, but rather part of a guided adaptation process that has many different solutions. We may have to consider this complexity to understand concentration adjustments in biological cells.

By using different objective functions for each target, it will be quite possible to combine different goals for targets in a multiple-output signaling system. Interestingly, a self-organizing cellular signaling system may also incorporate adaptive controls as objective functions, in particular when ion channels or membrane transporters are targeted. In this case, the goal of the system is to respond to an input (e.g. ISO at beta-2 receptors) such that the magnitude of another input such as calcium (ion channels as targets) is controlled. Signal transmission efficiency for such a target is not to be optimized to a maximum value, but instead to the appropriate ratio between the inputs. We have not followed up on this idea, but it may be worth further investigation as an example of natural computation for adaptive control.