## Abstract

A mathematical understanding of regulation, and, in particular, the role of feedback, has been central to the advance of the physical sciences and technology. In this article, the framework provided by systems biology is used to argue that the same can be true for molecular biology. In particular, and using basic modular methods of mathematical modelling which are standard in control theory, a set of dynamic models is developed for some illustrative cell signalling processes. These models, supported by recent experimental evidence, are used to argue that a control theoretical approach to the mechanisms of feedback in intracellular signalling is central to furthering our understanding of molecular communication. As a specific example, a MAPK (mitogen-activated protein kinase) signalling pathway is used to show how potential feedback mechanisms in the signalling process can be investigated in a simulated environment. Such ‘what if’ modelling/simulation studies have been an integral part of physical science research for many years. Using tools of control systems analysis, as embodied in the disciplines of systems biology, similar predictive modelling/simulation studies are now bearing fruit in cell signalling research.

- biochemical reaction network
- dynamics
- feedback
- pathway
- regulation
- systems biology

There is nothing more practical than a good theory.

David Hilbert

## Systems biology

As the name suggests, systems biology is a merger of systems theory, and molecular and cell biology. This conception of systems biology originates from Mesarović in 1968 [1], and it has received renewed interest following the advances that have been made with technologies to observe intra- and inter-cellular processes. The central dogma of systems biology is that it is the dynamic interactions of molecules and cells that give rise to biological function. The biological agenda of systems biology is defined by the following two questions relating intra- and inter-cellular processes within a cell and in cell populations:

How do the components within a cell interact, so as to bring about its structure and realize its functioning?

How do cells interact, so as to develop and maintain higher levels of organization and function?

The present article reviews the system-theoretical perspective on systems biology, focusing on the concept of feedback.

## Systems modelling

Systems biology, as defined above, considers signal- and systems-oriented approaches [2]. In this context, the present section first introduces the concepts of a system and signals, and argues the case for differential equation modelling. Systems theory tells us that in order to understand a dynamic system, that is the behaviour or functioning of a system, we need to perturb it systematically. Given a set of mathematical equations, in order to be able to identify parameter values from experimental data, the system has to be stimulated and a response recorded. A stimulus–response system is defined as a set of objects that are in relation to each other. More formally we write *S*⊆Ω×Γ, where the two objects Ω and Γ are sets of elements corresponding to some stimulus, and respectively the response to this stimulus. The symbol × denotes the set-product combination of different sets. A relationship is, in mathematical terms, a subset, denoted by ⊂, of a set product. The definition of a complex system now arises naturally as a system of systems. For a dynamic system, the elements of Ω and Γ are signals, e.g. observed concentration profiles

While it may be possible that different stimuli lead to the same response, we allow each stimulus to correspond to only one response profile. The mathematical relationship then becomes a mapping

The model (1) of a general dynamic system encompasses most, if not all, formal models that are considered in the engineering and physical sciences, applied mathematics, and molecular and cell biology. Whether we model with differential equations, automata, numerical values or symbols, and regardless of whether we choose a deterministic or stochastic framework, model (1) is the unifying abstraction. The aim for modelling at this level is the discovery of universal or generic organizing principles, that are independent of a particular realization, parameter values, cell line or organism.

While the stimulus-response model (1) is a suitable reflection of experimental reality (observing the temporal response to some stimulus), the aim of modelling is to hypothesize about its internal structure that lies between the input and output spaces of a system. An internal structure for the model is realized through the introduction of another object *X*, called the state-space of the system. Thus for (1), we have the following diagram expressing the relationships among objects of the state-space representation:

The principal aim of systems biology is to understand causal entailment in cells. Causation is the explanation of change, but, while changes occur in the realm of matter, causation is modelled as a between changes of states. For anything to change, either space or time, or both, have to be presupposed. A natural language to formalize the rate of change of system variables are differential equations:
where denotes the rate of change over time, which is often written as d*x*/d*t*, Θ is a vector of parameters (often fixed or constant over time), *x*(*t*) denotes the state vector of the system and *f* is a mapping that relates the change to the current state. As simple as eqn (2) appears, the abstraction that is behind this simplicity means that virtually all models of pathways that have been published, regardless of the number of molecules and their forms involved, are examples of this one equation. Take, for example, the very simple model of an enzyme kinetic reaction
where S refers to the substrate, E denotes the enzyme, SE is the substrate–enzyme complex and P is the product. Referring to the law of mass action, one arrives at the state-space model

The state vector of this system summarizes respectively the concentrations of the enzyme, substrate, substrate–enzyme complex and product *x*=(E, S, SE, P). The parameter vector is Θ=(*k*_{1}, *k*_{2}, *k*_{3}). Looking at the first equation in (3), it says that that the rate of change of the enzyme increases linearly with a decrease in SE and an increase in E. Further below, the enzyme kinetic reaction serves as a subsystem for signalling pathways. Such a modular approach is characteristic of control systems theory. From our brief introduction to systems theory, we can summarize the primary tasks in dynamic pathway modelling as:

Realization Theory: characterize model structures that could realize given stimulus–response data sets.

System Identification: determine values for model parameters, using experimental data or simulation studies.

Control Analysis: predict the consequence of changes to a pathway; in particular, modifications to parameters, cross-talk and the introduction/removal of feedback loops.

## Feedback regulation and control

Most dynamic motifs, i.e. oscillations, adaption, tracking and multi-stable switching behaviour, is the result of feedback interactions among system variables. The simplest model for the growth of molecular populations is for a constant stimulus, leading to a linear accumulation (integration). We can represent this by a signal-oriented block diagram:

where the response signal *y*(*t*) is mathematically formulated as

If the stimulus *u*(*t*) is a step-change, then *y*(*t*), the integral of *u*, is a straight line (Figure 1, top). The symbol of the block depicts this. Since the signal *y*(*t*) is unbounded, this model is rather unrealistic for any intracellular process. Cell functions are realized through feedback-regulated/controlled processes. We note that, in order to control, regulate, adapt, maintain or co-ordinate whatever, there must exist (implicitly) a goal or objective. To realize an objective, the current state of the system must be fed back to influence the change of state. In our simple ‘toy’ model of growth, we add a negative-feedback loop to the block diagram:

Note that in this block diagram, lines correspond to signals and the minus sign indicates the subtraction of the two values (e.g. concentrations). The equation for the response signal can be read off the block diagram as

As shown in Figure 1, the added negative-feedback loop has a dramatic effect on the dynamics, stabilizing the signal *y*(*t*).

## The effect of protein translocation

An important feature of feedback in biological systems is that it is often subject to a delay owing to the time taken to translocate molecules. This has important consequences for the stability of a process. As an example, here we consider the consequences of time delays in feedback loops in greater detail. In [3], an example is given for the JAK (Janus kinase)/STAT (signal transducer and activator of transcription) signal transduction pathway. In cell signalling, important mechanisms to transduce and relay signals are dimerization and (de-)phosphorylation. Upon stimulation of the membrane-bound JAK receptor, monomeric STAT5 is recruited towards the receptor and is then phosphorylated. This is followed by dimerization and eventually migration into the nucleus, where it binds to the promoter of target genes. This is followed by dephosphorylation and export back into the cytoplasm, where the molecules are returned to the process. This recycling out of the nucleus was generally considered to be a very slow process compared with the signal transduction into the nucleus. In [3], experimental data were initially fitted to a model that ignored the time delay occurring from the export back into the cytoplasm. The model generated smooth predictions, while the experimental data exhibited far more dynamic changes. The model was then extended to account for the time delay. The extended model revealed that STAT5 undergoes rapid nucleocytoplasmic cycles, continuously coupling receptor activation and target gene transcription, thereby forming a remote sensor between nucleus and receptor. The work in [3] is important, as it illustrates the crucial role of time delays and their relative size. An example for delays affecting dynamics in development can be found in [4].

To illustrate the effect of time delays in feedback loops, let us consider the simple systems with feedback considered above. We introduce in the feedback loop a time delay to simulate the type of effect discovered previously [3]. In the block diagram, *T*_{d} defines the translocation delay. These signal-oriented block diagrams are well-established in the engineering sciences [5].

As shown in Figure 1, the consequence of having a time delay, *T*_{d}, in the feedback path is to introduce damped oscillations in the system response. Thus time delays in feedback paths are generally destabilizing. The larger the time delay, the greater the destabilizing effect, in the limit leading to sustained oscillations.

## Dynamic modelling in cell signalling

The textbook idea of a signalling pathway is that of a linear cascade. Consider, for example, the model of a MAPK (mitogen-activated protein kinase) cascade [6]:

In the diagram, *x*_{1} corresponds to activated MKKK-P (phosphorylated MAPK kinase kinase), *x*_{2} to MKK-P (phosphorylated MAPK kinase), *x*_{3} to MKK-PP (double-phosphorylated MAPK kinase), *x*_{4} to MAPK-P (phosphorylated MAPK), *x*_{5} to MAPK-PP (double-phosphorylated MAPK). Inactivated forms are denoted with a tilde, ~. In signalling, activation/inactivation of proteins corresponds to phosphorylation/dephosphorylation, while, in some cases, one considers double-phosphorylations:

Denoting the kinase as *U*, the phosphatase as *P* and the protein as *X*, and assuming a constant phosphatase, a simple biochemical model of a signalling step is given by
where *v*_{1}(·) and *v*_{2}(·) are mappings, describing the reaction rates for phosphorylation and dephosphorylation respectively. We write for the non-phosphorylated form of the protein *X*, *u* for the kinase *U*, *p* for phosphatase *P* and *x* corresponds to the activated protein *X**. Referring to a power-law representation, one would have
where we choose a=b=c=1, such that

If we assume that the total =(*t*)+*x*(*t*) is constant for all time *t*, we require only this one differential equation to model a signalling step. See [7] for a more comprehensive study of this model of a signalling step. For the entire model of the generic MAPK pathway above, we now read off the reaction rates directly from the diagram:

The differential equations are then derived from the diagram as

In a full form we have for the model:

In addition, it is assumed that the following conservation relationships hold true:

The signal-oriented block diagram for a signalling step is readily obtained from the differential equations:

There is a new connector introduced here: the × indicates the product of two signals. There are also two blocks with parameters *k*_{1} and *k*_{2}, which are simple gains, for which the incoming signal is multiplied by these parameters. For more complex systems, the unambiguous encoding of a mathematical model with a signal-oriented block diagram allows us to collapse the previous diagram into a single block to represent the (de-)phosphorylation of through kinase *u*:

For clarity of representation, the constant will not be shown if the assumption of conservation law is clear. For double phosphorylation, let us consider the activation of *x*_{5} in terms of a signal-oriented block diagram:

This again we can collapse into a single block, without loss of information:

The linear MAPK pathway can then be represented by the compact block diagram:

This system is then another example of the state-space representation (eqn 2), where for *u*(*t*), we might assume an exponential decay *u*(*t*)=*e*^{−λt}, where an initial concentration of ligands is depleted through binding to the receptors on the cell surface. The greater the value of λ, the faster the ligands bind with receptors to form complexes. In [7], these series-connected subsystems of (de-)phosphorylation steps have been used to analyse pathways for their dynamic properties. In particular, the authors derived an expression for the signalling time, defined as the average time to activate a protein in the pathway, and the signal duration, characterized by an integral of the concentration profile.

In [6], a very similar model to the one above is modified by introducing a negative-feedback loop between the end product MAPK-PP and a Ras/MKKK complex at the top of the pathway. Kholodenko et al. [8] showed how ultrasensitivity, leading to switch-like behaviour, combined with negative feedback, can bring about sustained oscillations in this pathway. Approaches to test for feedback loops have been presented in [8,9]. In [10⇓–12], computational studies of feedback effects on signal dynamics in a detailed MAPK pathway model are presented.

## Feedback in signalling pathways

We now extend our discussion of feedback by introducing feedback loops from a protein *x*_{j} further down the pathway, up to *x*_{i}. We have two options illustrated by the following diagrams:

On the left-hand side, feedback manifests itself as a multiplicative term in the differential equation for *x*_{i}:
where for the function *F*(*x*_{j}), a commonly used choice is
where *n*≥1 defines the steepness of the feedback function, leading to ultrasensitivity of the (sub-)system. The subscript *I* of *K*_{I} stands for ‘inhibition’. The main requirement for the choice of a function *F*(*x*_{j}) is that at *x*_{j}=0, we should have *F*(*x*_{j})=1. Mechanistic interpretations and experimental evidence for these functions are discussed in [13⇓⇓–16]. Note the distinction between a mechanistic (or physical) and an operational (or phenomenological) definition for an interaction. An operational definition is based on observations, not necessarily requiring an interpretation/understanding of the physical interactions of the molecules involved, as would be the case for a mechanistic definition of kinetic behaviour [17]. For the feedback indicated on the right-hand side, and represented by *G*(*x*_{j}), there is an additional contribution to the activation of *X*:
leading to the following modified differential equation model

If *G*(*x*_{j}) is monotonically increasing with *x*_{j}, we have positive feedback, and, vice versa, if *G*(*x*_{j}) is monotonically decreasing in *x*_{j}, we are dealing with negative feedback. In a signal-oriented block diagram, we recognize the two situations as follows:

The above generalized motifs help in abstracting the system for the application of systems theory that has a rich history of half a century.

## Feedback regulation of a MAPK pathway

Based on the mathematical models introduced above, we here investigate feedback loops in the Ras/Raf/MEK [MAPK/ERK (extracellular-signal-regulated kinase)] kinase/ERK pathway [18,19] In the Ras/Raf/MEK/ERK pathway, corresponding to the generic MAPK pathway diagram in the previous section, Ras is the G-protein, Raf-1 is the MAPKKK, MEK is the MAPKK and ERK is the MAPK [16,20]. The feedback loop considered in the pathway is illustrated by the following schematic diagram:

The generic MAPK pathway can be recognized in the sequential activation of Raf-1, upstream near the cell membrane, followed by activation of the proteins MEK and ERK through structural modifications in the form of phosphorylations indicated by the Ps. The final double-phosphorylated ERK translocates into the nucleus of the cell, where it affects the transcription of genes. Double-phosphorylated ERK also phosphorylates RKIP and thereby releases Raf-1 from the Raf-1–RKIP complex, and the released Raf-1 in turn activates MEK. This positive-feedback loop can lead to various types of interesting pathway dynamics. For a mathematical model, let us consider a version with single phosphorylation steps as illustrated in Figure 2.

The variables of this system are *x*_{1}=Raf-1, *x*_{2}=MEK, *x*_{3}=ERK, *x*_{4}=RKIP. We first consider the pathway without any feedback loop. Although there may be grounds for a discussion whether the assumptions underlying a Michaelis–Menten model are appropriate in the context of signalling pathways, we here follow the literature [6,8,15] and choose the following set of equations:
where it is assumed that the conservation relationships _{1}=_{1}(*t*)+*x*_{1}(*t*), _{2}=_{2}(*t*)+*x*_{2}(*t*) and _{3}=_{3}(*t*)+*x*_{3}(*t*) hold true. Next, we consider the positive-feedback loop introduced by RKIP and which is denoted by *x*_{4}. The phosphorylation and dephosphorylation are described as before,
where _{4}=_{4}(*t*)+*x*_{4}(*t*). Note that *x*_{3}, activated ERK-PP, is acting on the phosphorylation of *x*_{4} (RKIP). The inhibitory effect of RKIP on the phosphorylation of *x*_{2} (MEK) is reflected by a change to the rate equation of *x*_{2}:
where *K*_{P} is a constant that defines the strength of the feedback and *p* defines the steepness of the response curve. The inhibitory effect, justifying the bar at the end of the arrow in Figure 2, is encapsulated in the square brackets. As the inactive form _{4}=_{4}−*x*_{4} increases, the value inside the square brackets decreases in the product with *x*_{1}(*t*). But because an increase in *x*_{3} leads to a decrease in _{4}, the overall feedback effect of *x*_{3} on *x*_{2} and the pathway is of a positive nature. The negative feedback from *x*_{3} (ERK-PP) to *x*_{1} (Raf-1) leads to an insertion in the equation for d*x*_{1}/d*t*:

Here *K*_{N} and *n* have the same role as *K*_{P} and *p* above, only that *n* and *N* stand for ‘negative’ rather than ‘positive’. For all proteins involved, conservation relationships hold for a constant total of the activated and non-phosphorylated form.

As shown in Figure 3, with only positive feedback added to the pathway and no transport delay, the pathway can display switch-like behaviour. The parameter values used in this simulation study can be found in Table 1. Switching dynamics have been found in various intracellular systems (e.g. [14,21]). While a negative-feedback loop was observed to have a stabilizing effect, a positive-feedback loop sharpens the response, making it ultrasensitive. Because the positive-feedback loop affects only proteins from MEK downwards, the Raf-1 concentration profile is not changed. Considering a negative-feedback loop, no transport delay and without the positive-feedback loop in the system, we observe that negative feedback can destabilize the response (see Figure 4). What can also be observed are lower steady-state values for Raf-1 and ERK.

Once a model is established, simulation allows quick studies of changes to the elements and parameters. For example, one way to make the model more realistic is to introduce a time delay between ERK near or inside the nucleus and its feedback effect on Raf-1 further up the pathway. In Figure 5, we introduce a transport delay in the negative-feedback loop with *T*_{d}=10 min. It can be observed that transport delays in a feedback loop (e.g. nucleo-cytoplasmic export) lead to increased oscillatory behaviour, turning the damped oscillations into sustained oscillations. Our next experiment is to change the feedback indices *n* and *p*, which were also introduced above and which define the sharpness or sensitivity of the feedback effect. In Figure 6, we find that, without transport delay, an increase from *n*=1 to *n*=2 in the negative-feedback loop also leads to sustained oscillations, but at a different frequency.

Our study demonstrates that there are various sources of sustained oscillations: negative feedback combined with ultrasensitivity [6], combined negative and positive feedback, as well as transport delays in negative-feedback loops. Oscillations have been investigated in various systems (e.g. [6,13,22]), and have been of interest in mathematical modelling for some time (e.g. [23,24]). An interesting question is to ask whether our model applies to a single cell or a population of cells. The role of feedback in intracellular dynamics has been investigated for some time in the literature (e.g. [4,25⇓⇓⇓–29]), and will, no doubt, play an important role in systems biology. Under the headings of phase-space analysis, bifurcation analysis, stimulus–response curves and stability analysis, systems theory provides mathematical and graphical tools to characterize the conditions under which certain dynamic motifs occur.

## Discussion and conclusions

To a large extent, the motivation for modern biomedical research is to elucidate the causal consequences of processes at the molecular level on the physiology and disease of organisms (Figure 7). At each level, feedback systems ensure a balanced response to environmental signals. This may be a regulatory response against perturbations or a controlled reaction to stimulus. The co-ordination of such multilayered and hierarchical processes remains an area of research with many interesting research problems.

The areas of genomics and bioinformatics have provided us with catalogues of components and their molecular characterization. Using technologies such as microarrays, two-dimensional gels and mass spectrometry, we are able to identify variables in intracellular processes and select them for modelling. It is for systems biology to describe and analyse the (dynamic) interactions of these variables (see [21,30⇓–32] for reviews and surveys). We presented the area of systems biology as a merger of systems theory with molecular and cell biology. The concepts of a system and signal were introduced with an application to signal-oriented block diagrams. These generalized motifs, which have a rich history in the engineering and physical sciences, serve as building blocks for systems of systems [33]. Albeit restricted to differential equation models, signal-oriented block diagrams are unambiguous representations compared with cartoon-like pathway maps. In an example, a simulation study revealed various possible sources for oscillations in signalling pathways. Transport delay, induced by the nucleo-cytoplasmic translocation of molecules, was highlighted as an important feature.

While there are already some examples that demonstrate the usefulness of mathematical modelling in the area of cell signalling, we should be under no illusion that it is possible to build accurate and precise physical models of molecular interactions. In cell signalling studies, experiments generate an indirect picture of processes within the cell. Quantitative immunoblotting, for example, involves various intermediate steps before a protein concentration can be quantified: (i) preparation of cellular lysate, (ii) immunoprecipitation, (iii) gel electrophoresis, (iv) transfer to a filter, (v) immunoblot, and (vi) imaging. Needless to say that the step from the *in vitro* to the *in vivo* situation provides a source for further speculation. The complexity of the systems under consideration, the technicalities of the experiments and the non-linear nature of interactions in biochemical reaction networks make it necessary to use mathematical modelling. In many interdisciplinary projects, it is the modelling process itself that is more important than the model. The discussion between the experimentalist and the theoretician, and the decision as to which variables to measure and why, form the basis for successful interdisciplinary collaborations. In the light of the interesting challenges this exciting area provides, systems biology is the art of making appropriate assumptions. The modelling process and the model are the means by which to complement the biologist's reasoning – no more, but no less either.

## Acknowledgments

M.U. and O.W. acknowledge the support of the U.K. Department for the Environment, Food and Rural Affairs (DEFRA), and the collaboration with the Veterinary Laboratories Agency (VLA), Weybridge, U.K. Work on the Ras/Raf/MEK/ERK pathway model is a joint effort with Walter Kolch from the Cancer Research U.K. Beatson Laboratories in Glasgow, Scotland. K.-H.C. acknowledges support by the Korean Ministry of Science and Technology through grants M10309000006-03B5000-00211 (Korean Systems Biology Research) and MG05-0204-3-0 (the 21C Frontier Microbial Genomics and Application Center Program). P.W.'s contribution was funded by Science Foundation Ireland under grant 03/RP1/I383. S.N.S. gratefully acknowledges the support of US NIH (National Institutes of Health) grant P20 CA 112963, US NIH grant K25 CA 113133, and Case Provost Opportunity Fund.

## Footnotes

Systems Biology: Will it Work?: Focused Meeting held at the University of Sheffield, U.K., 12–14 January 2005. Organized by M. Williamson and W. Blackstock (Sheffield). Edited by M. Williamson. Sponsored by BBSRC (Biotechnology and Biological Sciences Research Council), British Biophysical Society and GlaxoSmithKline.

**Abbreviations:**
ERK, extracellular-signal-regulated kinase;
JAK, Janus kinase;
MAPK, mitogen-activated protein kinase;
MEK, MAPK/ERK kinase;
MKK, MAPK kinase;
MKKK, MAPK kinase kinase;
-P, phosphorylated;
-PP, double-phosphorylated;
STAT, signal transducer and activator of transcription

- © 2005 The Biochemical Society