Research

An evolutionary approach for understanding of cellular response modulation

Manal A. Tashkandi 1, Mohammed Y. Refai 1, Lina A. Baz 2, and Mohammad Mobashir 3*

1 Department of Biochemistry, College of Science, University of Jeddah, Jeddah-21589 Saudi Arabia.

2 Department of Biochemistry, College of Science, King Abdulaziz University, Jeddah-21589 Saudi Arabia.

3 Department of Biomedical Laboratory Science, Faculty of Natural Sciences, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway.

* Correspondence: m.mobashir@cdslifesciences.com (M.M.)

Citation: Tashkandi, MA, Refai, MY, Baz, LA, and Mobashir, M. An evolutionary approach for understanding of cellular response modulation. Glob. Jour. Bas. Sci. 2025, 1(6). 1-10.

Received: January 07, 2025

Revised: February 28, 2025

Accepted: March 29, 2025

Published: April 21, 2025

doi: 10.63454/jbs20000029

ISSN: 3049-3315

Volume 1; Issue 6

 Download PDF File

Abstract: In T-cell signaling, the receptor, after receiving the signal activates downstream molecules that transmit the signal to the nucleus to control cellular functions such as proliferation, differentiation, and apoptosis. During this process the change in the protein-levels, post-translational modification of the proteins, and interaction strengths affect the final response of the cell. The objective is to develop mathematical model to understand the modulation of cellular response for signaling networks. We used signaling networks of fixed topology to summarize the collective understanding of flow of signaling information within a cell. To investigate the possible design principle, we developed an in-silico model using evolutionary algorithm and ordinary differential equations. To investigate the effect of interaction strength and concentration levels of the signaling molecules, the mutations were allowed in the protein-protein interactions and in the concentration levels of signaling molecules during the evolutionary period. From our results, we conclude that the kinetics of the final response of signaling pathway is predominantly controlled by the expression levels of proteins,  the variations in protein-levels within the cell plays critical roles in controlling the kinetic behavior of the cell and consequently the cellular functions, and an increase in the interaction strengths up to a certain level leads to strong activation patterns, i.e. in the cells become sensitive to stimuli. Thus, it appears that the kinetic properties of the signaling network mainly determine the response threshold of a cell while the quantitative behavior is controlled by the expression levels of the proteins.

Keywords: Signal transduction networks; modulation of cellular responses; fluctuations in protein concentration levels, kinetic parameter modification; mass-action kinetic model

1. Introduction

Signal transduction is a significant process for the cellular communication[1]. Based on the transduced signal strength from receptor to the effector signaling molecules through different signaling pathways, the cells undergo to different cellular processes like proliferations, apoptosis, differentiations, and divisions[2,3]. In signal transduction processes, an external stimulus is transformed into a cellular response through a network of proteins that produce different types of final response and accordingly alters the function and behavior of the cell[4-6]. Different forms of input-output relationships shown by biological signal transduction are known from experimental work. The change in the input signal strength, kinetic parameters, or the concentration levels of the signaling molecules can give rise to sustained, adapted responses, partially adapted, or oscillatory. Based on the response types the cells take the decision to specialized functions such as proliferation, differentiation, and apoptosis[7-9].

For better understanding about our theoretical model in terms of the biological pathway, we take MAPK (mitogen-activated protein kinase) pathway as an example, where TGFα binds to the EGFR (epidermal growth factor receptor), and then EGFR bound with TGFα activates the Grb2-Sos complex which activates the inactive Ras (bound with GDP) to active Ras (bound with GTP). This active Ras activates Raf, afterwards active Raf activates MEK and the active MEK activates the ERK. These activated ERKs translocate to the nucleus, where they phosphorylate and regulate various transcription factors leading to the changes in gene expression. This mitogen-activated protein kinase (MAPK) pathway shows different response types depending on the cell type and/or stimulus.

As we have seen in the previous paragraph that this MAPK pathway involves Raf, MEK (MAPK/ERK kinase), and ERK (extracellular signal-regulated kinase) and is considered to be centrally involved in cellular decision making processes where small quantitative differences often lead to major phenotypic changes. It has been shown that the upstream molecules induce quantitative and qualitative differences in the duration and magnitude of ERK activity that regulate the function and behavior of a cell. Recently, it has been shown that in case of peripheral T cells the ERK transient response leads to apoptosis and the sustained response to proliferation. But in other previous works, it has also been shown that depending on the different cell lines these response types have been correlated in different ways with different cellular functions such as proliferation, differentiation, and apoptosis. Here, we have given one example, if we interpret these conclusions in terms of cellular change by considering the previously shown data then we correlate it as very strong interactions will help the engineered PC12 cells to differentiate irrespective of the concentration levels fluctuations (sustained output response) and in the cellular systems with good interaction strengths and fluctuations in the concentration levels will block the differentiation in the wild-type PC12 cells (transient output response). The MAPK pathway is a prototype for the general scheme of signal transduction, in which, after receiving a signal from ligand-bound receptors, the involved proteins are altered (“activated”) by post-translational modifications. Subsequently, the active form activates other inactive proteins by means such as recruitment to specific locations, altering the enzymatic activity, or conformational changes exposing binding sites for further binding partners. To predict the function of a signaling module it is necessary to understand the design principles of signaling networks (SNs) that underlie the behavior and function. From experiments, neither the topology of a SN nor the kinetic parameters of its underlying elementary interactions are known in detail such that it remains open how sensitive the function of a network is to these parameters. And addition to it, it is also not possible to measure the concentration levels of all the mediator signaling molecules with in the cell. Therefore, it seems to be a challenging to explore the evolution of signaling networks allowing mutations of kinetic parameters and in the concentration level of the signaling molecules to mimic the subsequent acquisition of additional regulatory layers.

Recently, we have presented a theoretical approach where we investigated the role of kinetic parameters, addition of new nodes, and the input signal strength on the cellular response. In this work, we had already shown that these three factors do not play any role in producing transient cellular response but always produce sustained response. And additionally, it has also been shown that stronger interaction strengths and addition of new nodes produce strong and sustained response in the presence of input signals. In the absence of input signal there is no response at all. In the other previous studies, various modeling approaches have already been applied to investigate the behavior of SNs. Santos SDM et al., (2007)  has shown that the growth factor determines the topology of the MAPK signalling network and that the resulting dynamics govern the cell fate[10-13]. Francois and Hakim (2003) evolved genetic circuits to produce a variety of functional behaviors and demonstrated the vital role of post- transcriptional interactions, i.e. protein-protein interactions controlling gene regulation. This evolutionary approach has been extended by others to protein-protein interaction networks with specific functional characteristics: oscillators, bistable switches, homeostatic systems, and frequency filters. Yet, none of these approaches investigated in how far the formation of transient protein-protein complexes influences the generated networks or whether association, dissociation, or catalytic rates are critical for the network properties[14-21].

By using the same approach, we have investigated the SNs by keeping fixed but difference in the concentration levels of the receptor, mediator, and effector molecules and also by allowing the mutation in the concentration levels during the evolutionary phase. We found that the difference in the concentration levels of the signaling molecules play major role in controling the nature of the output response (sustained or transient) and the detailed parameter values are not critical for the functional response of the network. The interaction strength influences the sensitivity of the network, i.e. whether to respond or not to an input of given strength, rather than whether to respond with a transient or sustained output.

2. Results

2.1. Fundamental of biological pathway and evolution of RSNs: Signal Transduction event is mediated by the receptor, that transduces the extracellular signals by initiating a wide array of intracellular signaling pathways. One of the most studied signaling pathway is MAP kinase pathway where MAP kinase cascade communicates signals from active form of Ras i.e., RasGTP which phosphorylate Raf and Raf activates MEK and active MEK leads to the activation of ERK. Depending on the response type (transient or sustained) of fully active ERK the cell either go to proliferation or apoptosis or differentiation (Figures 1-4).  

Figure 1. The fitness (Fnorm) during the evolution of STNs for three different regimes of maximal kinetic parameters: weak (k < 10), moderate (k < 30), and strong (k < 100) interactions (Number of generations: 200, Number of STNs: 200, Threshold level: 1/10th of the initial concentration of protein). For all the systems, six different strength input signal were used which are 10^5, 10^3, 10^2, 1, 10, and 100. The concentration for the receptor, mediator, and effector molecules are 10, 5, and 1, respectively. (A) Evolution of STNs and (B) kinetics of the evolved STNs for three different kinetic parameter regimes in generation 1, 50, and 200. g1, g50, and g200 stand for generation 1, 50, and 200, respectively.

  

We evolve the SNs by assuming that the basic task of signal transduction is to provide an above-threshold response to a signal which is generated by a ligand binding to its receptor (see Methods) and the response is measured at a pre-selected node. When the activation state of this node crosses a threshold for an arbitrary time, the SN is considered to be successful in the detection of the corresponding signal and the fitness is increased by 1 () otherwise the fitness for the signal in question remains 0. For the task to detect a set of multiple signals of different strengths one by one the fitness contributions are summed up for all the input signals. In our simulations, signals are either present throughout the simulation or given as a pulse of fixed duration.

To investigate the modulation of cellular response from adapted to non-adapted and vice-versa, we evolved the networks by applying evolutionary algorithm under four different conditions. (i) Evolution of SNs due to variations in the kinetic parameters over generations while the concentration levels remain fixed. But the concentration levels of the receptor, mediator, and effector molecules are different (System I) (Figure 1A). (ii) Evolution of SNs due to variations in the concentration level of the receptor molecules and kinetic parameters over generations (System II) (Figure 2A). (iii) Evolution of SNs due to variations in the concentration level of the mediator signaling molecules and kinetic parameters over generations (System III) (Figure 3A). (iv) And evolution of SNs due to variations in the concentration level of the effector molecules and kinetic parameters over generations (System IV) (Figure 4A). Each system of System I has three different subsystems subsystem a (interaction strength < 10), b (interaction strength < 30), and c (interaction strength < 100). The other three systems (II, III, and IV ) also have three different subsystems as subsystem a (interaction strength and concentration level < 10), b (interaction strength and concentration level < 30), and c (interaction strength and concentration level < 100). The evolution with these four types of simulation conditions are investigated independently and compared using replicates with identical parameter settings but different seeds for the pseudo-random number generator. Then we measure the fitness (Fnorm) of the networks during the evolutionary period (Figure 1A, 2A, 3A, and 4A). Only the successfully evolved networks were allowed to enter next generation.

Figure 2 The fitness (Fnorm) during the evolution of STNs for three different regimes of maximal kinetic parameters: weak (k < 10), moderate (k < 30), and strong (k > 30) interactions (Number of generations: 200, Number of STNs: 200, Threshold level: 1/10th of the initial concentration of protein). For all the systems, 6 different strength input signal were used which are 10^5, 10^3, 10^2, 1, 10, and 100.

Here, we mainly looked over the effect of strong and weak interactions (System I) as well as the effect of the occurrences of change in the interaction strengths and in the concentration levels of the molecules involved in signal transduction during the evolutionary period (System II, III, and IV). In System Ia, when the interaction strength remains below a certain threshold (k < 10 in our model setting) the networks are unable to reach maximum fitness and at the same time the fitness of the population fluctuates significantly (Figure 1A). If the interaction strength remains below a threshold of 30 (System Ib) (Figure 1A), then the networks population reaches almost maximum fitness, but exhibits a considerable amount of fitness fluctuations (Figure 1A). Further increase in the interaction strength (k < 100) suppresses fluctuations in the fitness, i.e., a population has evolved in which virtually all the networks are able to detect every single input signal (System Ic) (Figure 1A). The evolution pattern of SNs in System IV where the mutation in the concentration levels of the effector molecules in addition to kinetics parameters has been allowed (Figure 4A) appears similar to the evolution pattern of the SNs in System I for all the three kinetic parameter regimes (Figure 1A). In Systems II and III, when the interaction strength remains below a certain threshold (k < 10) the networks are unable to reach maximum fitness (Figure 2A and 3A). In comparison to System Ia (Figure 1A) and System IV a (Figure 4A), the fitness of the population of Systems II and III fluctuates rapidly (Figure 2A and 3A). If the interaction strength remains below a threshold of 30, then the networks population reaches almost maximum fitness, but exhibits a rapid fitness fluctuation (System IIb and IIIb) (Figure 2A and 3A). Further increase in the interaction strength (k < 100) suppresses fluctuations in the fitness but still shows higher fluctuations than the systems Ic (Figure 1A) and IV c (Figure 4A). In comparison to systems I (Figure 1A) and IV (Figure 4A), systems II (Figure 2A) and III (Figure 3A) shows very high fluctuation in the fitness curve throughout the evolution period for all the three interaction strength regimes.  

2.2. Fixed but difference in concentration levels of the signaling molecules also acts as a possible factor for switching the sustained to transient cellular response: We have run the simulations for the evolution of SNs at constant concentration levels and with the possibility of variation in the interactions strengths in each generation for three different sets of concentration levels.  (a) Simulation with the concentration level set receptor molecule as 10, mediator molecule as 5, and effector molecule as 1 (System AI). (b) Simulation with the concentration level set receptor molecule as 5, mediator molecule as 10, and effector molecule as 1 (System BI). (c) Simulation with the concentration level set receptor molecule as 1, mediator molecule as 5, and effector molecule as 10 (System CI). 

Figure 3(A) Evolution of STNs and (B) kinetics of the evolved STNs for three different kinetic parameter regimes in generation 1, 50, and 200.

Each simulation were performed for three different regimes of interaction strength as we have mentioned in the previous section. We analyzed kinetics of the evolved SNs in System AI, BI, and CI (concentration levels are fixed but not equal) and observed the output node with transient response in all these systems for all the three different interaction strength regimes (Figure 1B). In this work, we have shown one of the result (System AI). This result appears completely different from our previously published result (where all types of signaling molecules have equal and fixed concentration level) where the output response is always  sustained. From this observation, we conclude that the difference in the concentration levels of the signaling molecules are responsible for producing the transient cellular response.

2.3. Effect of variation in the protein concentration levels on the kinetics of the cellular response: For the purpose to reveal more details about the effect of concentration level variations in the receptor, mediator, and effector molecules, we have also analyzed the kinetics of the other three systems (System II, III, and IV) and found that in these three systems too, fluctuations in the concentration level lead to transient response for all the three regimes of kinetic parameters with few exceptions in case of very strong interactions (Figure 2B and 3B). In case of very strong interactions most of the networks show transient response but there are also considerable number of networks which switch their transient response to sustained response. In all the other three systems (II, III, and IV), we observed final output response as transient (for most of the evolved SNs). During the middle of the evolutionary period, final output response of most of the evolved networks show transient nature but after full evolution there are also considerable number of networks with switch like and strongly activated sustained response (Figure 2B, 3B, and 4B).

In our previous work, we have already shown that in the presence of equal and fixed concentration levels, the networks do not show any kind of transient behavior when the interaction strength varies during the evolution while in systems I, II, III, and IV, after 10 generations, we observed the transient activation. From this observation we can conclude that the strength of input signals and the protein- protein interaction strengths do not play any role either together or independently in controlling the transient nature of the final output response but stronger input signals and stronger interaction strengths lead to the strong and sustained activation pattern. The possible factors for the transient response adaptation are variations in the concentration level of the signaling molecules and difference between the concentration levels of receptor, mediator, and effector molecules. 

Figure 4 Evolution of STNs by mutating kinetic parameters and the concentration of the effector molecule A3. (A) Evolution of STNs and (B) kinetics of the evolved STNs for three different kinetic parameter regimes in generation 1, 50, and 200.

2.4. Effect of variation in the protein concentration levels on the kinetics of the cellular response: The activation patterns evolving when strong interactions are permitted appear independent of the input signal strength, a behavior also known for the MAPK cascade in certain systems. We extended our analysis to investigate the dose-response of the evolved networks beyond the range of signaling strength used for selection. Before evolution for all the four simulating systems, at very weak input signal, the networks do not produce a considerable response with respect to the threshold (Figure 5A, B, C, and D). After a few generations, all the networks are able to evolve and show strong activation even at very low input signal strengths (Figure 5A, B, C, and D). Yet, the networks require a signal to become active confirming that the evolution did not generate self-activating networks, which would not be excluded by our choice of the fitness function. Thus, the strength of the input signal does not affect the activation pattern of the SN. Therefore, it appears  that the generic behavior of a SN is switch-like when facing the task to ’somehow’ detect a signal. For all the systems (Figure 5), 10 different strength input signal were used which are 10−9, 10−7, 10−5, 10−3, 10−2, 0.1, 1, 10, 50, and 100. We have also analyzed the partially active nodes which always show transient activation and evolved the SNs where the input signals remain present for a defined time period only but we get the results as shown in our previous work.

Figure 5 Dose-response relationship of the evolved STNs for moderate kinetic parameter for all the Systems The graph shows the maximum response of the output node of the best network versus the input signal in the respective generations 1, 100, and 200. For all the systems, 10 different strength input signal were used which are 10^9, 10^7, 10^5, 10^3, 10^2, 0:1, 1, 10, 50, and 100.

3. Discussion

In this work, we have discussed four different simulation conditions and investigated the evolution of SNs where the primary task of signal transduction is to detect a signal without pre-determining a desired kinetic. In these simulations, mutations were allowed during the evolutionary process in the kinetic parameters and the concentration levels of the receptor, mediator, and the effector molecules. The simulations, where the mutations were allowed only in kinetic parameters and the concentration levels of the receptor, mediator, and effector molecules are fixed but not equal in amount to each other (System I), has the generic solution with a transient activity of the output node as long as the signal is present. Irrespective of the interaction strength almost all the networks are able to detect the signal and produce the final response as transient response. While the networks evolved at weak interaction strength also have a considerable number of networks which respond differently and only to those signals which are sufficiently strong (in strength) and the response does not cross the threshold level for signals with very weak strength, in brief we conclude that the network respond only to some of the signals (Figure 1B). In a cellular system with weak interactions could therefore probably not provide reliable cellular decisions. And the other three systems where the mutations were allowed also in the concentration levels of the molecules involved in signaling also show the generic solution with transient response with few exceptions. Similar to System I, these three systems (System II, III, and IV) respond in a similar way in terms of the response kinetics (transient/partially adapted response).

The evolved networks in a stationary population show stable activation pattern against strong interactions and the transient response in System I, suggesting that this type of response is the generic cellular behavior when the presence of a signal is sufficient information for a cell. We observe that in System I the response kinetics does not alter after about 30 − 50 generations and maintains transient nature but the kinetic parameters still change. While in case of also change in the concentration levels (System II, III, and IV), the majority of the SNs show transient response within 100 generations but afterwards the response kinetics fluctuate too much and a large number of SNs change their transient behavior to sustained behavior but the kinetic parameters and concentration level still change. We interpret these outcomes in the following way: when the networks population with enhanced fitness appear they give rise to multiple clones that have largely similar kinetic parameters which is similar to the bottleneck effect, i.e. many networks do not generate offsprings due to their low fitness and only similar networks pass on to the next generation. Following this phase the networks start to diversify again and a large range of the allowed kinetic parameter regime is explored. The diversification also indicates that there is a large number of solutions in the parameter space to ’solve’ the fitness function. These solutions appear to be connected in a large set as the different kinetic parameters can be explored by the SNs without losing their fitness. Alternatively, one can view this situation as overfitting as the quite large number of parameters allows the networks to ’solve’ the fitness function in many different ways.

We also compared the first system (System I) (where we have simulated the SNs in the presence of six different input signal strengths and during the evolutionary process mutations in the kinetic parameters (interaction strengths) were allowed but the concentration level of signaling molecules were kept fixed (but the concentration levels of effector, mediator, and effector molecules are not equal) throughout the evolutionary period) with our previous work (here, we have simulated the SNs in the presence of six different input signal strengths and during the evolutionary process mutations in the kinetic parameters (interaction strengths) were allowed but the concentration level of signaling molecules were kept fixed (but the concentration levels of effector, mediator, and effector molecules are equal) throughout the evolutionary period). In all the evolved networks in System I, we found that irrespective of the interaction strengths the nature of final response is always transient (Figure 1B). While in our previous work (where the concentration level of signaling molecules were kept fixed (and the concentration levels of effector, mediator, and effector molecules were also equal) the strength of output response is directly proportional to the interaction strength and nature of response is always sustained. It means the final response of the evolved networks becomes stronger and stronger until they reach to the maximum and almost all the networks show switch like response. From this comparative analysis we conclude that the difference in the concentration levels of the receptor, mediator, and effector molecules is responsible for this transient response (Figure 1B, 2B, 3B, and 4B).

As we have shown in our previous work, the variation in the kinetic parameters but with fixed and equal concentration level of the signaling molecules shows the generic cellular response as sustained[14,15]. So, in conclusion from the current and previous works, from System I, we conclude that the difference between the concentration levels of the receptor, mediator, and effector molecules produces the transient response although we do not allow mutation in the concentration levels during the evolution phase. Or in another way we can say that irrespective of the protein-protein interaction strengths, difference between the concentration levels of signaling molecules (receptor, intermediate, and effector) acts as a possible factor for controling the transient response when the mere detection of the presence of a signal is relevant (System I, II, III, and IV). Previous studies demonstrated that the variation of kinetic parameters and addition of nodes is sufficient to evolve the networks that have a defined output response[23,28-32]. However, this required corresponding fitness functions that encode the mathematical property of the desired system, e.g., bistability. These models can only be applied to those systems which are known to have such behavior, but often the exact behavior of the SNs is not known. Therefore, the creation of a fitness function that encodes the task that a cell solves under certain experimental conditions, may be more beneficial in determining possible and likely behavior of the underlying SNs.

Our current model evolves SNs with the variations in interaction strength and concentration levels of the signaling molecules of the SNs. The kinetic parameters and the concentration levels of all the signaling molecules are in general hard to access experimentally by the experimentalists. In the simulations, we see that the exact value of the kinetic parameter plays a minor role. Unfortunately, the number of possible parameter sets to generate a certain behavior is large such that the topology alone is not likely to predict the function of the network reliably. However, it is not the only interaction strengths and input signal which vary but also protein concentration levels affect the behavior of a SN and protein concentrations bear the advantage of being experimentally quite easily accessible variables. Furthermore, stimulation of receptors virtually never occurs in isolation and therefore the interaction of signal transduction pathways can become relevant. The questions of how does cross-talk affect the network’s behavior and how does it affect the evolution of the SN is worth pursuing using our approach by embedding a SN into the wider context of co-evolving signaling networks [34-36].

As we saw in the previous subsection that the strength of input signals and the protein-protein interaction strengths do not play any role either together or independently in controlling the transient nature of the final output response but interaction strengths leads to the strong and sustained activation pattern but from the results it appears that the concentration level variation as one of the factor leading output response to the transient behavior. Since the interaction strength and input signal strengths do not have any direct impact on the nature of the final response but when we evolved the networks allowing the variations in the protein concentration levels during evolution then we found that there are majority of the evolved networks which show transient response. So we can say from these observations that the fluctuation in protein concentration levels definitely have some kind of roles in controlling the transient output response.

4. Conclusions
From these observations, we reach to the conclusions: first point, irrespective of the protein-protein interaction strengths, the variation in the concentration levels signaling molecules (receptor, intermediate, and effector) produces transient response when the mere detection of the presence of a signal is relevant (System II, III, and IV). Second, from System I, we conclude that the difference between the con- centration levels of the receptor, mediator, and effector molecules also produces the transient response although we do not allow mutation in the concentration levels during the evolution phase. So from these two points we can say in summary as mainly the difference between the concentration levels are responsible for the transient responses.

We used our recently published model to represent a signal transduction pathway allowing two post- translational modifications of similar to the MAPK cascade. In order to transduce the signal, we have included protein-protein interactions, protein phosphorylation and dephosphorylation. Double phosphorylated proteins act as fully activated and single phosphorylated molecules as partially activated molecules. Note, that the term phosphorylation is used for convenience as any other post-translation modification adding a small chemical group, lipid, protein or carbohydrate modifying a protein’s spectrum of interaction partners or enzymatic activity are covered by the model. The interaction between the signaling proteins are set up randomly to create the initial population which were used to evolve under the conditions of variation in either only kinetic parameters or variations in both kinetic parameters and the concentration level of the signaling molecules. Initially, all the proteins are inactive. One of the initially present inactive signaling molecules is designated as receptor and receives the signal to become activated. Once the receptor receives the signal then it can activate other signaling molecules. All the signaling molecules are allowed to phosphorylate or dephosphorylate each other and the final products will be formed depending on the complex. All the reactions in this model are bimolecular, autophosphorylation and homodimer formation are not allowed. Every molecule that becomes partially (single phosphorylated) or full active (double phosphorylated) can interact with any other molecule in any state. These interactions lead to complex formation. The complexes can dissociate without changes to its constituents or upon modifying on of it by means of phosphorylation or dephosphorylation. The interaction of two partially active molecules produces either one of them being fully activated (dual phosphorylated) or inactivated (dephosphorylated) without changing the other reacting partner’s state (attributing it an enzymatic role). Which of the possible reactions are realized is determined randomly once at the beginning (with constraints, see next paragraph), thus setting up the network topology. In addition, one randomly selected double phosphorylated protein, different from the receptor node, is designated the output node (effector molecule). It represents the molecule which will eventually induce the cellular response.

5. Methods

5.1. Systems set up descriptions: In the result section, we have shown the data mainly from four systems and each of these systems have three different subsystems for three different ranges of kinetic parameters (system 1Aa, 2a, 3a, and 4a for k < 10, system 1Ab, 2b, 3b, and 4b for k < 30, and system 1Ac, 2c, 3c, and 4c for k < 100). In addition to these four systems we have two more systems for the fixed concentration levels as systems 1B and 1C. The systems 1A, 1B, and 1C have fixed concentrations levels throughout the evolutionary period have the concentration level sets for signaling (receptor, mediator, effector) molecules as (10, 5, 1), (10, 5, 1), and (10, 5, 1), respectively [34-36].

5.2. Dynamical mass action kinetics equation: A network (Figure 6) consists of the above mentioned signal transduction pathway where input signal (green), inactive (blue), partially active proteins (cyan), fully active (red) proteins and the protein complexes (yellow) are represented as nodes. The interaction matrix (Aij) between all these molecules including complexes are represented as +1 (production/generation), −1 (degradation/dissociation), and 0 (no interaction). The entries of Aij are chosen once for a network under the constraints that the total amount of each protein is conserved and the SN generated does have a stable inactive state. The entries Aij generate the index/indices for the reactant(s) xp(r) for the reaction r. Each arc encoded by the interaction matrix is associated with a weight that represents the kinetic parameters with which production or degradation takes place. The dynamics xi of the node i is governed by the equation:

equation.pdf    ……….(1)

here, rtot is the total number of reactions and kr denotes the kinetic parameter of the reaction number r. And k as dimensionless in the sense that the time is scaled appropriately and the concentrations are normalized such that the numeric values of first- and second-order reactions approach a similar range. 

The fitness (Fnorm) of a SN was tested by calculating its response to ns = 6 different signals. For every signal n, the dynamics of the pre-selected output node is tested whether it exceeds a threshold f. This threshold  is defined to be the relative fraction of the dual phosphorylated protein to the total 10 amount of the protein. If the output node crosses the threshold f at any point during the dynamics the network gains a fitness contribution Ffactor(n) = 1. The normalized fitness Fnorm is calculated as the average fitness contribution for all signals which are weighted equally:

equation_1.pdf equation_2.pdf

Hence, the maximum fitness of a network will be 1 when it detects all signals or 0 when there is no above-threshold response to any of the signals. Different strength input signals used for calculating the kinetics are 10−5, 10−3, 10−2, 1, 10, and 100. But for Dose-Response relation curve, we have used ten different strength input signals as 10−5, 10−4, 10−3, 10−2, 0.1, 1, 5, 10, 50, and 100.

5.3. Evolutionary algorithm: We applied an evolutionary algorithm to optimize and evolve signaling networks under defined constraints. The approach was designed to explore the parameter space of protein interactions while maintaining a fixed network topology. Before initiating the evolutionary process, we generated a diverse set of networks to serve as the initial population. Specifically, we created N = 200 networks, each constructed using the same randomly generated interaction matrix. This matrix defined the connectivity among three proteins, with each protein modeled in three distinct phosphorylation states: unphosphorylated, monophosphorylated, and dual-phosphorylated. Although the interaction matrix was identical across all networks, diversity was introduced by assigning different randomly selected kinetic parameters to each network, as described in Equation 1.

The evolutionary process was governed by two mechanisms of variation. First, mutation of the kinetic parameters was allowed, whereby reaction rate constants (kr) were perturbed to explore alternative dynamic behaviors. Second, variations in protein concentration levels were permitted, enabling the networks to adapt through changes in molecular abundance. Importantly, the algorithm did not allow the addition of new nodes, thereby preserving the structural integrity of the original interaction matrix throughout the evolutionary process【15,33,34】. For the pre-evolved signaling networks, the kinetic parameters were sampled randomly from a uniform distribution within the range 0.001 to 1.0, ensuring sufficient variability in the initial population.

For each network in the population, a normalized fitness measure (Fnorm) was computed to evaluate its performance relative to the desired signaling behavior. Based on these fitness values, elite selection was performed, whereby the most successful networks were retained for reproduction. The selected networks were then subjected to mutation, either through changes in kinetic parameters or through variations in protein concentration levels, as outlined above. To maintain population size, each successful network was replicated into four copies, thereby populating the subsequent generation while keeping the total number of networks constant at N = 200.

This evolutionary cycle was repeated for 200 generations, allowing gradual refinement of network properties and convergence toward optimized signaling behaviors. At each generation, the population size remained fixed, ensuring controlled evolutionary pressure and comparability across iterations. The system dynamics underlying each network were modeled using ordinary differential equations (ODEs) derived from Equation 1. These ODEs were numerically solved using MATLAB 7.9.0Author Contributions: Conceptualisation, M.A.T., M.Y.R., L.A.B., and M.M.; software, M.M.; investigation, M.M.;  writing—original draft preparation, M.A.T., M.Y.R., L.A.B., and M.M.; writing—review and editing, M.A.T., M.Y.R., L.A.B., and M.M.; visualisation, M.A.T., M.Y.R., L.A.B., and M.M.; supervision, M.M.; project administration, M.M. The author has read and agreed to the published version of the manuscript.

Funding: Not applicable.

Acknowledgments: We are grateful to the Department of Biochemistry, College of Science, University of Jeddah, Jeddah-21589 Saudi Arabia and Department of Biomedical Laboratory Science, Faculty of Natural Sciences, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway for providing us all the facilities to carry out the entire work.

Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.

Data Availability Statement: All the related data are supplied in this work or have been referenced properly.

References

  1. B.N. Kholodenko, Cell-signalling dynamics in time and space, Nature Reviews Molecular Cell Biology. 7 (2006) 165–176. doi:10.1038/nrm1838.
  2. Q. Cui, Y. Ma, M. Jaramillo, H. Bari, A. Awan, S. Yang, et al., A map of human cancer signaling, Molecular Systems Biology. 3 (2007). doi:10.1038/msb4100200.
  3. Z. Culig, F.R. Santer, Androgen receptor signaling in prostate cancer, Cancer Metastasis Rev. 33 (2014) 413–427. doi:10.1007/s10555-013-9474-0.
  4. J.E. Purvis, G. Lahav, Encoding and Decoding Cellular Information through Signaling Dynamics, Cell. 152 (2013) 945–956. doi:10.1016/j.cell.2013.02.005.
  5. X. Xia, M.S. Owen, R.E.C. Lee, S. Gaudet, Cell-to-cell variability in cell death: can systems biology help us make sense of it all? 5 (2014) e1261–11. doi:10.1038/cddis.2014.199.
  6. B.E. Housden, N. Perrimon, Spatial and temporal organization of signaling pathways, Trends in Biochemical Sciences. 39 (2014) 457–464. doi:10.1016/j.tibs.2014.07.008.
  7. M.K. Morris, J. Saez-Rodriguez, P.K. Sorger, D.A. Lauffenburger, Logic-Based Models for the Analysis of Cell Signaling Networks, Biochemistry. 49 (2010) 3216–3224. doi:10.1021/bi902202q.
  8. J. Saez-Rodriguez, L.G. Alexopoulos, M. Zhang, M.K. Morris, D.A. Lauffenburger, P.K. Sorger, Comparing Signaling Networks between Normal and Transformed Hepatocytes Using Discrete Logical Models, Cancer Research. 71 (2011) 5400–5411. doi:10.1158/0008-5472.CAN-10-4453.
  9. J. Saez-Rodriguez, L.G. Alexopoulos, G. Stolovitzky, Setting the Standards for Signal Transduction Research, Science Signaling. 4 (2011) pe10–pe10. doi:10.1126/scisignal.2001844.
  10. 10. S.L. Spencer, S. Gaudet, J.G. Albeck, J.M. Burke, P.K. Sorger, Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis, Nature. 459 (2009) 428–432. doi:10.1038/nature08012.
  11. 11. J. Rameseder, K. Krismer, Y. Dayma, T. Ehrenberger, M.K. Hwang, E.M. Airoldi, et al., A Multivariate Computational Method to Analyze High-Content RNAi Screening Data, Journal of Biomolecular Screening. 20 (2015) 985–997. doi:10.1177/1087057115583037.
  12. 12. K. Moelling, K. Schad, M. Bosse, S. Zimmermann, M. Schweneker, Regulation of Raf-Akt Cross-talk, J. Biol. Chem. 277 (2002) 31099–31106. doi:10.1074/jbc.M111974200.
  13. 13. W. Kolch, M. Halasz, M. Granovskaya, B.N. Kholodenko, The dynamic control of signal transduction networks in cancer cells, Nature Reviews Cancer. 15 (2015) 515–527. doi:10.1038/nrc3983.
  14. 14. M. Mobashir, T. Madhusudhan, B. Isermann, T. Beyer, B. Schraven, Negative Interactions and Feedback Regulations Are Required for Transient Cellular Response, Sci. Rep. 4 (2014). doi:10.1038/srep03718.
  15. 15. M. Mobashir, B. Schraven, T. Beyer, Simulated evolution of signal transduction networks, PLoS ONE. 7 (2012) e50905.
  16. 16. U.S. Bhalla and Ravi Iyengar, Emergent Properties of Networks of Biological Signaling Pathways, Science. 283 (1999) 381–387. doi:10.1126/science.283.5400.381.
  17. 17. J.A. Papin, B.Ø. Palsson, Topological analysis of mass-balanced signaling networks: a framework to obtain network properties including crosstalk, Journal of Theoretical Biology. 227 (2004) 283–297. doi:10.1016/j.jtbi.2003.11.016.
  18. 18. L. Bardwell, X. Zou, Q. Nie, N.L. Komarova, Mathematical Models of Specificity in Cell Signaling, Biophysical Journal. 92 (2007) 3425–3441. doi:10.1529/biophysj.106.090084.
  19. 19. K.P. Hofmann, C.M.T. Spahn, R. Heinrich, U. Heinemann, Building functional modules from molecular interactions, Trends in Biochemical Sciences. 31 (2006) 497–508. doi:10.1016/j.tibs.2006.07.006.
  20. 20. H.A. Harrington, E. Feliu, C. Wiuf, M.P.H. Stumpf, Cellular Compartments Cause Multistability and Allow Cells to ProcessMore Information, Biophysical Journal. 104 (2013) 1824–1831. doi:10.1016/j.bpj.2013.02.028.
  21. 21. W.W. Chen, M. Niepel, P.K. Sorger, Classic and contemporary approaches to modeling biochemical reactions, Genes & Development. 24 (2010) 1861–1875. doi:10.1101/gad.1945410.
  22. 22. J.J. Hornberg, B. Binder, F.J. Bruggeman, B. Schoeberl, R. Heinrich, H.V. Westerhoff, Control of MAPK signalling: from complexity to what really matters, Oncogene. 24 (2005) 5533–5542. doi:10.1038/sj.onc.1208817.
  23. 23. L. Qiao, R.B. Nachbar, I.G. Kevrekidis, S.Y. Shvartsman, Bistability and Oscillations in the Huang-Ferrell Model of MAPK Signaling, PLoS Comput Biol. 3 (2007) e184. doi:10.1371/journal.pcbi.0030184.st002.
  24. 24. D.Y. Wang, L. Cardelli, A. Phillips, N. Piterman, J. Fisher, Computational modeling of the EGFR network elucidates control mechanisms regulating signal dynamics, BMC Systems Biology. 3 (2009) 118. doi:10.1186/1752-0509-3-118.
  25. 25. C. Sun, R. Bernards, Feedback and redundancy in receptor tyrosine kinase signaling: relevance to cancer therapies, Trends in Biochemical Sciences. 39 (2014) 465–474. doi:10.1016/j.tibs.2014.08.010.
  26. 26. P.J. Roberts, C.J. Der, Targeting the Raf-MEK-ERK mitogen-activated protein kinase cascade for the treatment of cancer, Oncogene. 26 (2007) 3291–3310. doi:10.1038/sj.onc.1210422.
  27. 27. J.D. Thomas, T. Lee, N.P. Suh, A FUNCTION-BASED FRAMEWORK FOR UNDERSTANDING BIOLOGICAL SYSTEMS, Annu. Rev. Biophys. Biomol. Struct. 33 (2004) 75–93. doi:10.1146/annurev.biophys.33.110502.132654.
  28. 28. S.D.M. Santos, P.J. Verveer, P.I.H. Bastiaens, Growth factor-induced MAPK network topology shapes Erk response determining PC-12 cell fate, Nature Cell Biology. 9 (2007) 324–330. doi:10.1038/ncb1543.
  29. 29. R.J. Orton, M.E. Adriaens, A. Gormand, O.E. Sturm, W. Kolch, D.R. Gilbert, Computational modelling of cancerous mutations in the EGFR/ERK signalling pathway, BMC Systems Biology. 3 (2009) 100. doi:10.1186/1752-0509-3-100.
  30. 30. E. Nikonova, M.A. Tsyganov, W. Kolch, D. Fey, B.N. Kholodenko, Control of the G-protein cascade dynamics by GDP dissociation inhibitors, Mol. BioSyst. 9 (2013) 2454. doi:10.1039/c3mb70152b.
  31. 31. M.R. Birtwistle, J. Rauch, A. Kiyatkin, E. Aksamitiene, M. Dobrzyński, J.B. Hoek, et al., Emergence of bimodal cell population responses from the interplay between analog single-cell signaling and protein expression noise, BMC Systems Biology. 6 (2012) 109. doi:10.1214/aos/1176346577.
  32. 32. T. Tian, A. Harding, How MAP kinase modules function as robust, yet adaptable, circuits, Cell Cycle. 13 (2014) 2379–2390. doi:10.4161/cc.29349.
  33. 33. K. Kaneko, Shaping robust system through evolution, Chaos. 18 (2008) 026112. doi:10.1063/1.2912458.
  34. 34. Almowallad S, Jeet R, and Mobashir M.   Systems-level understanding of toxicology and cardiovascular system. Glob. Jour. Bas. Sci. 2024, 1(1). 1-16.
  35. 35. S. Almowallad, R. Jeet, and M. Mobashir.  A systems pharmacology approach for targeted study of potential inflammatory pathways and their genes in atherosclerosis. Glob. Jour. Bas. Sci. 2024, 1(2). 1-12.
  36. 36. M. Mobashir, B. Schraven, T. Beyer (2012) Simulated Evolution of Signal Transduction Networks. PLoS ONE 7(12): e50905.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of Global Journal of Basic Science and/or the editor(s). Global Journal of Basic Science and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. 

Copyright: © 2025 by the authors. Submitted for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Loading

Bibliography