Self-organization of signal transduction

We propose a model of parameter learning for signal transduction, where the objective function is defined by signal transmission efficiency. We apply this to learn kinetic rates as a form of evolutionary learning, and look for parameters which satisfy the objective. This is a novel approach compared to the usual technique of adjusting parameters only on the basis of experimental data. The resulting model is self-organizing, i.e. perturbations in protein concentrations or changes in extracellular signaling will automatically lead to adaptation. We systematically perturb protein concentrations and observe the response of the system. We find compensatory or co-regulation of protein expression levels. In a novel experiment, we alter the distribution of extracellular signaling, and observe adaptation based on optimizing signal transmission. We also discuss the relationship between signaling with and without transients. Signaling by transients may involve maximization of signal transmission efficiency for the peak response, but a minimization in steady-state responses. With an appropriate objective function, this can also be achieved by concentration adjustment. Self-organizing systems may be predictive of unwanted drug interference effects, since they aim to mimic complex cellular adaptation in a unified way.


Introduction
Signal transduction systems are often modeled as networks of biochemical kinetic equations implemented as continuous-time dynamical models using differential equations [3,12].If we regard a subset of species as inputs, and make sure that the system always converges to equilibrium values by using weakly reversible equations [8,1,24], we may transform these models into a set of matrices fulfilling the role of input-output transfer functions, i.e. a mapping from sustained input signal levels to steady-state concentrations for all target species [18].Protein signaling functions (psfs) are a systemic generalization of individual dose-response functions, which are usually described by Hill equations [2].In contrast to Hill equations, which are not available for enzymatic reactions, which only calculate relative concentrations, and which only work for one reaction in isolation, the psf system calculates enzymatic and complex formation reactions in a complex systemic environment using absolute concentrations [18].In addition, the reaction times to equilibrium are calculated as delay values, and the dynamic shape ('transients') is also available for further analysis.
In this paper, we want to ask the question of optimality of signal transduction.From an evolutionary standpoint, we assume that any biological signal transduction system is constructed with optimized efficiency of signal transmission.Furthermore, we assume that cells have the ability to adapt to perturbations of protein concentrations and changes in extracellular signaling by reinstating signal transmission efficiency.In the following, we investigate this question using a biologically realistic system -beta-adrenergic signaling in a submembrane compartment of a mouse embryonic fibroblast-for a single input scenario, focusing on a selected target species as relevant output or actuator of the system (Fig. 1 A).
Experimental analysis of signal transduction systems has shown fold-change responses to changes in input [10,6].Accordingly, input-output transfer functions usually follow the shape of hyperbolic (saturating) curves, which are equivalent to sigmoids for logarithmically scaled input [18].In Fig. 1 B we show the effect of a knock-out (KO) for a RGS protein in an experimental assay in yeast [26], and compare this with the effect in the model system [4].We see that the effects of the RGS KO on dose-response signaling efficiency are robust across very different cellular systems.
Usually, when we use a computational model to investigate perturbations, we only study the effects as reflected in the simulation.By utilizing optimization in terms of signaling efficiency, we can make the system itself adjust to the perturbation.In this way, we are studying signal transduction as a self-organizing system, which uses objective functions to adapt.This basic idea is extremely powerful, and could be used with different kinds of constraints on parameters, reaction times, etc. and with different, multiple input scenarios for larger systems.To explore this question further is of significant importance in assessing cellular health and functioning.

Example System: GPCR signaling in a submembrane compartment
Fig. 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, [4]).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 [4].The dose-response transfer functions were derived as in [18].

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.1: Kinetic rates of the sample system, adapted from Sabio-Rk [25] and Brenda [17], adjusted to experimental cAMP time-series data [4] 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 (Fig. 1).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 steadystate 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 inputoutput 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 Fig. 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 (Fig. 1A).The optimization method used is a simplex algorithm [23].Fig. 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 ) in a biological system and under optimization define an optimal operating range for a signaling system and be of relevance in disease processes [18].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 (Fig. 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 and a wide variability in R C,S above a certain threshold in d is robust against different types of signaling input (cf.also Fig. 4A).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.[21,14]).As expected, this low variation did not significantly affect the quality of a set of kinetic rates (supplemental table).
We further analysed the parameter combinations with different signal transmission efficiency.In Fig. 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 [11].Standard deviation for all parameters in our case is similar (Fig. 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 selforganization 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.Fig. 3A), 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 downregulations, to maintain sensitivity of signaling, but reducing the receptor beta-2 leads to up-regulation.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 [22].Another idea would be to rank reaction systems defined by kinetic parameters as in Fig. 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 reuptake, 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 (Fig. 3B).As a result, a number of internal concentrations will change.Fig. 3C 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 [16], but also in cancer where intercellular signaling (e.g. by cytokines [15]) 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 steadystate (cf.2.2, Fig. 4A).This requirement of invariance for the steady-state is sometimes considered a form of 'robustness' or 'homeostatic regulation' of the cellular response [20,13], 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 (Fig. 4A).In Fig. 4B 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 Fig. 4C.The simplicity of Fig. 4C 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 selforganizing 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 to further investigate as an example of natural computation for adaptive control.

Discussion
The idea to look for parametric optimization as the basis of realistic cellular properties has been applied with success to metabolic fluxes [9].In that case, optimization of growth is usually regarded as the single objective function.Here we use another objective function, signal transmission efficiency, to study the adaptive response of a signaling system to perturbation in concentrations.This equals maximization of concentration change in a target species in response to input signal.Optimization of growth and optimization of signal transmission are therefore related.
Signal transmission in a cellular system may have multiple functions: transcription factor activation, which may be related to the cell cycle, to morphological change or to adjustment of protein concentrations, cytosolic kinase/phosphatase activation with multiple cellular targets, membrane protein activation such as ion channels or receptors, etc.We assume that signal transduction has been optimized by evolution, and that concentration adjustment exists to maintain effective signal transmission.We have shown how to optimize a single-input single-output system for both speed and signal efficiency due to the basic properties of kinetic equations [18].It is easy to extend the present discussion to optimize for multiple outputs in parallel, and a system could also be optimized for a number of I/O functions.In that case other measures, such as mutual information, may also be employed as objectives.We used optimization in a two-step process: (a) in evolutionary learning, in order to find kinetic rates for estimated concentrations (b) in cellular adaptation, in order to re-calibrate the model in response to perturbations in concentrations, changes in extracellular signaling, or to effect a transient vs. plateau-like response.This also means that we have transitioned from a biologically defined system that is built bottom-up from available parameters and data to a self-organizing, learning system that adjusts to changes on evolutionary or individual time-scales.
Matching models to experimental time-series data is an under-constrained problem that often yields large ranges of suitable parameters [19,11].Analysing signals and targets to find optimal transmission parameters may be used to further constrain and investigate the parametric space.Signal transmission efficiency in a biological system can be measured directly [7,5] and be compared to what is theoretically possible given a set of equations.During evolution, new protein subtypes develop with a different set of interactions and regulations.This corresponds to an adjustment of the available set of reactions, a type of structural learning to overcome the bottlenecks that are a result of tightly specified molecular kinetics and limited adjustment of concentrations.

Figure 1 :
Figure 1: A. Biochemical Reaction System with Selected Input(red) and Output(blue) B. (Top) ISO-pVASP transfer function (with RGS KO,red)) in the MEF model [4] (Bottom) Experimental Response to RGS KO for GPCR signal-response in a yeast model [26] C. (Top) Distribution of Extracellular Signals Used to Calculate Signal Transmission.From a baseline signal S 0 [10nM-900nM] the extracellular signaling level S t rises to 110%,150% and 200%, but not above 1µM.(Bottom) Transfer Functions for ISO → pVASP.Shown is the original model, optimization by delay, by R C,S , and by both.

Figure 2 :
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. Table

Table 2 :
Initial concentrations (in nM) of species in the sample system

Table 3 :
Reaction times (d S ) and response efficiency (R C,S