Targeted reduction of complex models with time scale hierarchy – a case study

Mathematical Biosciences 185 (2003) 123–152 www.elsevier.com/locate/mbs

Targeted reduction of complex models with time scale hierarchy – a case study Leiv Øyehaug a

a,1

, Erik Plahte

a,*

, Stig W. Omholt

b

Department of Mathematical Sciences, Agricultural University of Norway, P.O. Box 5035, 1432 A
s, Norway b Department of Animal Science, P.O. Box 5025, Agricultural University of Norway, 1432 A
s, Norway Received 28 January 2002; received in revised form 17 March 2003; accepted 20 June 2003

Abstract With the increasing flow of biological data there is a growing demand for mathematical tools whereby essential aspects of complex causal dynamic models can be captured and detected by simpler mathematical models without sacrificing too much of the realism provided by the original ones. Given the presence of a time scale hierarchy, singular perturbation techniques represent an elegant method for making such minimised mathematical representations. Any reduction of a complex model by singular perturbation methods is a targeted reduction by the fact that one has to pick certain mechanisms, processes or aspects thought to be essential in a given explanatory context. Here we illustrate how such a targeted reduction of a complex model of melanogenesis in mammals recently developed by the authors provides a way to improve the understanding of how the melanogenic system may behave in a switch-like manner between production of the two major types of melanins. The reduced model is shown by numerical means to be in good quantitative agreement with the original model. Furthermore, it is shown how the reduced model discloses hidden robustness features of the full model, and how the making of a reduced model represents an efficient analytical sensitivity analysis. In addition to yielding new insights concerning the melanogenic system, the paper provides an illustration of a protocol that could be followed to make validated simplifications of complex biological models possessing time scale hierarchies.
2003 Elsevier Inc. All rights reserved. Keywords: Melanogenesis; Time scale hierarchy; Singular perturbations; Bistability; Melanogenic switching

*

Corresponding author. Tel.: +47-64 94 8865; fax: +47-64 94 8810. E-mail address: [email protected] (E. Plahte). 1 Present address: Norwegian Defence Research Establishment, 2007 Kjeller, Norway. 0025-5564/$ - see front matter
2003 Elsevier Inc. All rights reserved. doi:10.1016/S0025-5564(03)00095-6

124

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

1. Introduction In the post-genomic era, biological research will become almost synonymous with the efforts to understand the functional expression of genes within the context of integrated biological systems. Biology is finally in a position to start revealing the causal links between genotype and phenotype in the wide sense. Non-linear dynamic models will play a key role in this process by which known feedback loop structures, signal transduction pathways and associated metabolic pathways will be combined into coherent models spanning several systemic levels. Considering the bewildering complexity displayed by biological systems, many of these models will have a very large number of state variables and a very complex interaction structure. Such models are normally not suitable for detecting generic regulatory principles, modular behaviour and essential causal features, especially when several models are linked together to form even larger models. Thus there is a growing demand for mathematical tools by which essential features of these models can be captured and formulated in simpler mathematical representations without loosing too much of the realism of the original models. Given the presence of a time scale hierarchy, application of singular perturbation techniques represents an existing elegant method of developing minimised mathematical representations that agree quantitatively with more complex models of the system. Application of singular perturbation methods in mathematical biology is particularly widespread in models of enzyme processes. In fact, singular perturbations are implicitly applied in all models applying Michaelis–Menten kinetics, see e.g. [1, Chapter 5]. These matters have been extensively studied for relatively simple enzyme reaction networks, e.g. [2–5] as well as for more complex enzyme reaction networks, e.g. [6,7]. The ideas originating in enzyme kinetics were suggested by Liao and Lightfoot [8] and Delgado and Liao [9] to work even for general metabolic pathways and some authors have applied these to concrete biological systems: Maki and Keizer [10] made a model for oscillatory insulin secretion analyseable by taking into account the time scale hierarchy. Similarly, assumptions on time scales enabled Tang et al. [11] to arrive at a tractable, analyseable model of calcium channel kinetics. In this paper we adopt the model reducing approaches of Maki and Keizer and Tang et al., but we will to a larger extent than these authors focus on the approximation properties of the simplified model solutions compared to the original model solutions. We (i) conduct a rigorous demonstration of how a large, complex model of a concrete biological system can be approximated by a substantially simplified model by singular perturbation methods without sacrificing biological realism and (ii) perform a detailed analysis of the resulting simplified model in order to examine the behaviour of the biological system in question. This is done by using a recently developed model of melanogenesis. Melanogenesis is the main material basis of the colour patterns observed in mammals and proceeds within pigment producing cells called melanocytes. Through the melanocyte stimulating hormone receptor (MC1-R) at the surface of the melanocyte, the melanocyte responds to variation in the concentration of external agents, especially melanocyte stimulating hormone (a-MSH) and its antagonist agouti protein, by producing varying amounts of the black/brown eumelanin or the yellow/reddish pheomelanin [12–15]. The two types of melanins are synthesised within organelles called melanosomes [16]. The melanosomes serve as Ôcolour packagesÕ which are transferred to surrounding keratinocytes. The basic regulatory mechanisms of this process have been the subject of intense study over the past three decades. Thus, a comprehensive amount of data concerning the various steps and

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

C1

DOPA

C2

DQ E1

DOPA DQ

E6

DHI E7 DHI Ox. IndQ

DC E3

GSHRB

Tyrosinase

Tyrosinase

Tyrosine

TRP2

DHICA

E2 E4 TRP1

LDC

E8

IndQCA E5

Eumelanin

SW1 GSHR GSH GSSG P1 SW2 Glu-DOPA P2

DQ DOPA

125

γ -GTP

CDOPA DQ P3 DOPA

CDQ P4 Pheomelanin

Fig. 1. The main components of the melanogenic system within individual melanosomes and interactions between them. The conversions labeled by C1 and C2 (see Eqs. (C1) and (C2) in Appendix A.1) constitute the initial common pathway of melanogenesis, the subprocesses labeled E1,. . .,E8 (Eqs. (E1)–(E8) in Appendix A.1) constitute the eumelanin pathway and the subprocesses labeled P1,. . .,P4 (Eqs. (P1)–(P4) in Appendix A.1) constitute the pheomelanin pathway. The processes labeled SW1, SW2 (Eqs. (SW1) and (SW2) in A.1) are crucial for the switching of melanin production between eumelanin and pheomelanin. Abbreviations are explained in Table 1.

actors in the eumelanin and pheomelanin pathways have accumulated (the most important ones are depicted in Fig. 1 and the corresponding chemical equations are given in Appendix A.1). In a recent paper [17] we tried to translate the available empirical knowledge into the language of nonlinear system dynamics. Comprising 18 coupled differential equations (see Appendix A.2 and A.3), this model is too complicated for many purposes. A reduction of a complex model like this one by singular perturbation methods will have to be a targeted reduction as one first has to decide which mechanisms, processes and features one wants to accentuate or disclose by the simplification. Normally, one wants to get a deeper understanding of the essential regulatory aspects of the complex model. However, what counts as essential aspects has of course to be defined in advance according to a set of explanatory goals. In melanin synthesis, an important problem is to explain why and how a kind of switch seems to operate when the melanocytes produce mixed melanins that are either dominated by pheomelanin or by eumelanin. In this perspective, the essential aspects will be the dynamics of the main regulating agents, the regulating mechanism, and the subprocess(es) responsible for this switch between what seems to be two clearly distinguished states of the melanocyte. Within the melanosome, eumelanin and pheomelanin are formed as the final products of the eumelanin and pheomelanin pathways, respectively (Fig. 1). The intermediate DOPAquinone (DQ) participates in the initiation of both these pathways, suggesting that DQ plays a key role in the switching of the melanocyte between producing mainly eumelanin and mainly pheomelanin. The hypothesis of Ito [16], recently supported by us [17], seems to be the most plausible candidate for a switching mechanism. Ito hypothesised that DQ binds itself covalently to the enzyme glutathione reductase (GSHR) (process SW1 in Fig. 1 and Eq. (SW1) in Appendix A.1). This enzyme catalyses the formation of glutathione (GSH) (process SW2 and Eq. (SW2)) which is the first metabolite in the pheomelanin pathway. Thus, at high levels of DQ, pheomelanin production is strangulated, whereas eumelanin production proceeds at a high rate. Conversely, when DQ levels are low, GSH reductase is not inactivated, and the pheomelanin pathway is initiated, leading to production of pheomelanin.

126

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

The hypothesis of Ito suggests that the dynamics of GSH reductase are essential in melanogenesis. We expect that the pattern formation process (i.e. the rate at which patterning is regulated by external agents) is very slow compared to most of the cellular subprocesses, like enzymatic and redox reactions leading to the synthesis and distribution of melanin. This motivates our hypothesis that there are three time scales in melanogenesis, namely 1. The time scale of the fast internal processes, i.e. the biochemical processes within the melanocyte. 2. The time scale of dissociation of the DQ GSH reductase complex (GSHRB in Fig. 1). 3. The time scale of external processes, i.e. the time scale of variation in the concentrations of agouti and a-MSH. Since the external processes are not an integrated part of the melanogenic system, the levels of external agents act as input to the intracellular system and since there is no feedback from within the cell to the external agents in the model, the time scales on which we focus represent intracellular processes. Mathematical models of systems comprising a time scale hierarchy typically express the presence of different time scales with one or several small parameters. Developing such models involves choosing a certain scaling (see Appendix B.1). The choice of scaling may be directed by what one knows about the reaction rates of the biological system in question, or by hypotheses on the system. For example, if there are strong indications that a certain subprocess is of particular importance for the behaviour of the system, we should choose a scaling according to the time scale of this subprocess. Since the process towards finding a scaled model involves targeting a specific subprocess, we will denote the scaled model a targeted mathematical representation (TMR
) of the biological system. Here,
is a measure of the difference in time scales in the biological system. When
is small, approximating
by 0 would simplify the model considerably. The resulting model will be referred to as TMR0 . It is essential that the TMR
is chosen so that the solution of the corresponding TMR0 is in fact an approximation to the solution of the TMR
for small values of
. This is the case provided that a number of assumptions on the wellbehavedness of the functions involved are satisfied (see Appendix B.2). Furthermore we should validate the TMR0 towards the TMR
as well as towards the real biological system. This validation should provide answers to the following questions: Does the TMR0 make sense biologically? Does it confirm our hypothesis on the essential subprocesses? Does it approximate TMR
in relevant parameter and variable domains? Does it give new insight into the original model and the real system, e.g. by disclosing robustness properties, generic properties and suggestions for future experiments?

2. Targeted model reduction of the melanogenic system 2.1. A possible switching mechanism The choice of scaling will be motivated by a conjecture on the switching mechanism. The switching could be caused by a number of different mechanisms, like e.g. steep sigmoidal interactions. But preliminary investigations lead us to suspect that the switching is a consequence of

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

127

bistability and a hysteresis effect due to changes in the a-MSH and agouti levels. Bistability requires a positive loop [18–20], and the most promising candidate seemed to be þ







GT !Gs !Dq !GT ; where GT represents the level of not inactivated GSH reductase, Gs is the level of GSH, Dq represents the DQ concentration and the signs and the arrows represent positive resp. negative effects. We expect this process to represent a key regulatory mechanism. 2.2. Choice of time scale In order to obtain a mathematical representation that targets the specific process which we suspect represents the key mechanism of the biological system, we must choose a scaling accordingly. Thus, our time scale of choice is t ¼ k
26 t, where k
26 is the rate of dissociation of the DQ GSH reductase complex, which we assume is small compared to most conversion rates. Using this scaled time we derive a model of melanogenesis which is of the form (B.2) (Appendix B.1). This choice of scaling is by no means the only possible, nor the only one worth considering to arrive at a TMR
. In the following the reader is advised to confer Table 2 for the definition of dimensionless parameters occurring in the model. The parameter
0 which occurs in the model of melanogenesis (see Eq. (2) below) is the ratio of the tyrosinase decay rate to the rate of dissociation of inactivated GSH reductase, i.e.
0 ¼ c=k
26 . Thus
0 is the ratio of a relative rate rA corresponding to a process A to a relative rate rB corresponding to a process B, i.e.
0 ¼ rA =rB . If process A is much slower than process B, then
0 is small. Table 1 The main substances involved in melanogenesis and the variables associated with them on dimensional and dimensionfree form Chemical component

Variable

Non-dimensional variable

a-MSH Agouti protein MSH receptor Agouti receptor Tyrosinase DOPA DOPAquinone (DQ) LeucoDOPAchrome (LDC) DOPAchrome (DC) DHICA IndQCA DHI IndQ Eumelanin Glutathione (GSH) Not inactivated GSH reductase Glu-DOPA CysteinylDOPA (CDOPA) CysteinylDOPAquinone (CDQ) Pheomelanin

Rm Ra M A T Dp Dq Lc Dc Da Ia Dh I E Gs GT Gp Cp Cq P

rm ¼ j1 Rm =k
1 ra ¼ j2 Ra =k
2 M A xT ¼ max sðMÞ=c xDP ¼ j19 Dp =k7 xDQ ¼ j19 Dq =k7 xLDC ¼ j8 Lc =k7 xDC ¼ j19 Dq =k7 xDA ¼ j19 Da =k7 xIA ¼ j19 Ia =k7 xDH ¼ j19 Dh =k7 xI ¼ j19 I=k7 e ¼ j19 E=k7 xGS ¼ j19 Gs =k7 xGT ¼ GT =EGSHR xGP ¼ j19 Gp =k7 xCP ¼ j19 Cp =k7 xCQ ¼ j19 Cq =k7 p ¼ j19 P =k7

128

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

Table 2 The non-dimensional parameters given in terms of the parameters of the biochemical and mathematical models Parameter name

Definition



0 l a rm ðt Þ ra ðt Þ b1 ðMÞ b2 ðMÞ b3 ðMÞ a c a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 h1 h2 h3 h4 ds de dp

k
26 =k7 c=k
26 k
1 =k
26 k
2 =k
26 j1 Rm ðk
26 tÞ j2 Ra ðk
26 tÞ sðMÞ= max sðMÞ hj3 ðMÞST j19 max sðMÞ=ðck72 Þ /j5 ðMÞ max sðMÞ=ðck7 Þ k~26 =k
26 k~25 j19 EGSHR =ðk72 Þ k13 j19 ETRP 1 =ðk72 Þ j22 =j19 j8 =j19 k11 j19 ETRP 2 =ðk72 Þ k15 =k7 k14 =k7 k17 j19 EDO =ðk72 Þ k18 =k7 j19 EGSHR =k7 k21 j19 Ec-GTP ðk72 Þ k23 =k p7ffiffiffiffiffiffiffi j19 Kqa =k7 j19 Kc =k7 j19 Kh =k7 j19 Kp =k7 cs =k7 ce =k7 cp =k7

The parameter
is defined as Relative rate of dissociation of inactivated GSH reductase k
26
¼ : ¼ Relative rate of conversion of DQ to leucoDOPAchrome k7 Whether
 1 or not remains an open question. It is known that the relative decay rate of tyrosinase is very small since the half life of this substance is approximately 10 h (see e.g. [16]). Thus the relative decay rate of tyrosinase is ln 2=36 000 s
1 . The relative rate of conversion of DQ to LeucoDOPAchrome is known to be 7.6 s
1 [16]. Hence

0  0:0000025. Based on this minute value we hypothesise that both
and
0 are small numbers;
;
0  1. 2.3. Development of a TMR
from the full melanogenic model In the model there is no feedback interaction from within the melanocyte to the cell surface. Thus we may consider the dynamics of MSH and agouti receptor at the surface of the melanocyte and the intracellular dynamics separately.

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

129

2.3.1. Dynamics at the surface of the melanocyte Introducing the time dependent (and dimension-free) variables rm and ra (defined in Table 1) for the concentrations of a-MSH and agouti in the surroundings of the cell and letting M and A represent the fractional occupancy of MSH and agouti receptors, respectively, Eq. (A.1) (Appendix A.2) are transformed to the non-autonomous, dimension-free equations dM ¼ rm ðt Þð1
M
AÞ
lM; dt dA ¼ ra ðt Þð1
M
AÞ
aA; dt

ð1Þ

where l and a are the ratios between the relative dissociation rates of MSH receptors and agouti receptors, respectively, and the relative dissociation rate of inactive GSH reductase k
26 . It has been argued that the number of melanocortin receptors expressed by the melanocyte is small (< 100 according to Robbins et al. [15]) rendering a continuous treatment of the receptor dynamics erroneous. On the other hand the estimate of Lu et al. [13] is much larger (103 –5 · 104 ). In this paper and in [17] we assume that the number of receptors is sufficiently large to defend a continuous treatment of the receptor dynamics. This is a convenient, but by no means critical, assumption, since the main focus of this paper is on the internal dynamics, not on the MSHreceptor dynamics. We presume that there is a positive, monotonic relationship between the fractional occupancy of activated MSH receptor (M) and the activity and synthesis of tyrosinase. On dimensional form the regulation of tyrosinase levels are described by the first of Eq. (A.2). In terms of dimensionfree variables and parameters the regulation is described by dxT ¼
0 ðb1 ðMÞ
xT Þ; dt

ð2Þ

where xT represents the concentration of tyrosinase, b1 ðMÞ represents the synthesis rate of tyrosinase, which we assume depends on M, and
0 was defined in the previous subsection. Eqs. (1) and (2) are decoupled from the rest of the model, i.e. the solutions of these equations merely act as input into the model of the intracellular model. The model we refer to as the TMR
therefore only involves the ODEs to be presented below, not Eqs. (1) and (2) above. 2.3.2. Melanogenic pathways within the melanosome The two initial steps of intramelanosomal melanogenesis are common for the eumelanogenicand the pheomelanogenic pathways, and they involve respectively the tyrosinase mediated hydroxylation of tyrosine to 3,4-dihydroxyphenylalanine (DOPA), and oxidation of DOPA to DOPAquinone (DQ) (subprocesses C1 and C2 in Fig. 1 and Eqs. (C1) and (C2) in Appendix A.1). The eumelanin pathway proceeds by cyclisation of DOPAquinone (DQ) to leucoDOPAchrome (LDC). A cascade of intermediate metabolites are subsequently formed and converted (subprocesses E1–E8 in Fig. 1 and Eqs. (E1)–(E8) in Appendix A.1) until eumelanin is formed as the final product. The pheomelanin pathway is initiated as DQ and GSH react to form Glu-DOPA, which through a sequence of reactions is converted to pheomelanin (processes P1–P4 in Fig. 1 and Eqs. (P1)–(P4) in Appendix A.1).

130

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

In addition, DQ may inactivate GSH reductase leading to a decrease in the GSH synthesis rate (subprocesses SW1 and SW2 in Fig. 1 and Eqs. (SW1) and (SW2) in Appendix A.1). We suspect that the inactivation of GSH reductase is the key regulator of melanogenic switching. Using the scaled time and dimensionless variables the time rate of change of the variables associated with the concentration of the main substances involved in melanogenic switching and the variables representing the concentrations of the metabolites of the common pathway and the eumelanin and pheomelanin pathways are described by dxGT ¼ 1
ð1 þ aIðxGT ; xDQ ÞÞxGT ; ð3aÞ dt dy ¼ gðt ; xGT ; y;
Þ; ð3bÞ dt where xGT represents the ratio between not inactivated GSH reductase and the total amount of GSH reductase in free and bound form, the variable xDQ represents the level of DQ and the function I is described below. The parameter a is the ratio between the rate of binding of DQ to GSH reductase and dissociation of inactivated GSH reductase and
was explained in Section 2.2. Furthermore, the vector y comprises the variables corresponding to the concentrations of the substances occurring in intramelanosomal melanogenesis and is defined in Eqs. (A.5), (A.7), (A.9) in Appendix A.3. Similarly, the right hand side g incorporates the interactions between the metabolites of intramelanosomal melanogenesis. This function is defined in Eqs. (A.6), (A.8), (A.10) in Appendix A.3. Eqs. (3a) and (3b) constitute our TMR
. The inactivating function I : ½0; 1Þ  ½0; 1Þ ! ½0; 1Þ measures the relative rate of GSH reductase inactivation as a function of xGT and xDQ . In [17] we tested the fit of the model to data using the functional relation xnDQ ; IðxGT ; xDQ Þ ¼ n xDQ þ dHn


where n and d are positive numbers that were estimated by a least-squares method and H could be either a constant or equal to xGT . Since we want our model to be flexible and capable of handling different functional forms for the modelling of GSH reductase inactivation, the function I is not specified at this stage, but, motivated by the results in our previous paper we assume that I has the following properties: 1. I is continuous and differentiable with respect to both arguments. 2. The speed of inactivation of GSH reductase increases with the level of DQ, i.e. I is increasing with xDQ . 3. The speed of inactivation of GSH reductase decreases with the level of available GSH reductase, i.e. I is decreasing with xGT . 4. There is no inactivation of GSH reductase when the level of DQ is zero, i.e. IðxGT ; 0Þ ¼ 0 for all values of xGT . Due to the complexity of the TMR
(3a) and (3b) an analysis of this model is not feasible. The next step in the process is therefore to apply singular perturbation methods to obtain an approximate model representing the melanogenic system in a realistic manner.

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

131

Proving that the solution of the full system (B.2), corresponding to the TMR
, approximates the solution of the reduced system (B.3), corresponding to the TMR0 , is not straightforward. There exist a number of results on the convergence of the reduced system solution to the solution of the full system when
! 0, all of them stating that the convergence is uniform granted that the full and reduced systems fulfill a number of assumptions on the wellbehavedness of the functions involved. For our purpose the Levin and Levinson theorem (see Appendix B.2) is used to prove that the TMR
can be approximated by the TMR0 . To the TMR
(3a) and (3b) corresponds the reduced system du ¼ 1
ð1 þ aIðu; vÞÞu; ð4aÞ dt ð4bÞ 0 ¼ gðt ; u; y0 ; 0Þ; where u and v denote the singular perturbation approximations of xGT and xDQ , respectively, and y0 is the singular perturbation approximation of y occurring in Eq. (3b), i.e. Eqs. (4a) and (4b) are obtained simply by setting
¼ 0 in Eq. (3b) and renaming the variables. We prove in Appendix C that the solution of the reduced system (4a) and (4b) is a uniformly valid approximation to the TMR
(3a) and (3b), i.e. we prove the following result: Result. Provided that
is small and the conjecture in Appendix C.1 holds, the solution of the reduced system (4a), (4b) is a uniformly valid approximation to the solution of the full system (3a) and (3b) in any closed and finite time interval where Assumption 2 (see Appendix B.2) holds. The Conjecture referred to in the above result states that all eigenvalues of the Jacobian of the rapid subsystem (Eq. (3b)) have negative real part whereas Assumption 2 in Appendix B.2 asserts that the right hand sides of the full system are continuous and differentiable. 2.4. Biological interpretation of the TMR0 As a consequence of the Result, one ODE (Eq. (4a)) and a number of easy-to-solve algebraic equations (Eq. (4b)) are sufficient to describe the intracellular dynamics of melanogenesis. The variable v, occurring in Eq. (4a), represents the level of DQ. It can be expressed in terms of u by solving Eq. (4b). This is straightforward, it suffices to add the two expressions in Eq. (A.6), setting the result equal to zero and thereafter solving gpð1Þ ¼ 0, where gpð1Þ is defined in Eq. (A.10). This gives the relation 1 vðu; bÞ ¼ ½ðb
ds
cuÞ þ ððb
ds
cuÞ2 þ 4ds bÞ1=2 ; ð5Þ 2 where b ¼ b2 ðMÞxT and the remaining variables and parameters are explained below. Thus, the TMR0 consists of Eqs. (4a) and (5). There are three parameters in the expression for v which means that there are four parameters in the ODE (4a) plus the parameters of the inactivating function I. All the variables and parameters in Eq. (4a) can be interpreted biologically. • The variable u, 0 6 u 6 1, is the ratio of the number of not inactivated GSH reductase molecules (in free or complex bound form) to the total number (not inactivated as well as inactivated) of GSH reductase molecules, which is conserved.

132

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

• The parameter a measures the ratio of the binding rate of DQ and GSH reductase to the dissociation rate of inactivated GSH reductase. Its value has been estimated to 25–30 [17]. • The quantity b, which appears as a constant in Eq. (5), is in fact the quantity b2 ðMÞxT . Here the variables M and xT are the, possibly time dependent, solutions of Eqs. (1) and (2) and represent the levels of activated MSH receptor and the activity of tyrosinase, respectively. Tyrosinase activity is considered to be the main controlling factor of melanin synthesis, emphasising the importance of b. In the analysis of the TMR0 (Section 3) we investigate the significance of b for the dynamical behaviour of the TMR0 , i.e. we study the impact of varying tyrosinase activities on the behaviour of the melanogenic system. • The product cu represents the rate of synthesis of GSH, i.e. c is the maximum rate of GSH synthesis. • Finally, the parameter ds represents the relative decay rate of GSH. The (non-dimensional) levels of eumelanin (e0 ) and pheomelanin (p0 ) depend crucially on the levels of GSH reductase and DQ; e0 ¼

1 vðu; bÞ; de

p0 ¼

c vðu; bÞ u ; dp vðu; bÞ þ ds

ð6Þ

where de is the relative decay rate of the non-dimensional eumelanin level e0 and dp is the relative decay rate of the non-dimensional pheomelanin level p0 . The expressions in Eq. (6) are obtained by solving Eq. (4b). If most of the GSH reductase is inactivated, i.e. u is small, the variable v is large and much eumelanin is synthesised. If the opposite is the case, i.e. little or no enzyme is inactivated, the value of u is close to 1, v is small and Eq. (6) indicates that little eumelanin is synthesised. On the other hand, provided that ds is small compared to v, the levels of pheomelanin are high.

3. Analysis of the TMR0 The TMR0 (Eqs. (4a) and (5)) is simple enough to be treated analytically. We are particularly interested in investigating the potential of the TMR0 to exhibit complex behaviour such as bistability. 3.1. Confined sets From the interpretation of u as a fraction of a conserved quantity u is required to be confined to the interval [0,1]. Inserting the endpoints of this interval into Eq. (4a) it is evident that uðt Þ stays within [0,1] when it is within the interval initially. A smaller confined set exists. Define the function g as gðu; bÞ ¼

1 1 þ aIðu; vðu; bÞÞ

and let g
ðbÞ ¼ minu2½0;1 gðu; bÞ and gþ ðbÞ ¼ maxu2½0;1 gðu; bÞ. By insertion it is observed that the interval bounded by g
ðbÞ below and gþ ðbÞ above is a confined set for Eq. (4a).

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

133

3.2. Critical points Critical points are determined as the solution(s) of the equation Kðu; bÞ ¼ 1
ð1 þ aIðu; vðu; bÞÞÞu ¼ 0:

ð7Þ

Critical points are stable if oK=ou is negative and unstable if oK=ou is positive. For a fixed value of b, Kð0; bÞ ¼ 1 > 0 and Kð1; bÞ ¼
aIð1; vð1; bÞÞ 6 0, hence the continuity of K ensures that there exists at least one critical point. If there is exactly one critical point this point must be stable since K then must be decreasing with u over the whole interval [0,1]. Unless there are points of multiplicity higher than one, there are an odd number of critical points (see e.g. Fig. 2(c) for an illustration). Numbering these points in increasing order, i.e. u1 < u2 <    < u2pþ1 , points of odd subscripts are stable and points of even subscripts are unstable. When there are 2p þ 1 distinct critical points, p þ 1 of these are stable, i.e. there are p þ 1 attractors in the system. In the case of three critical points there are two attractors and one repelling

(a)

(b)

(c)

1

1

1

0.5

0.5

0.5

0

0

0

−0.5 0

0.5

1

−0.5 0

0.5

−0.5 0

0.5

1

(e)

(d) 1

1

0.5

0.5

0

0

−0.5 0

1

0.5

1

−0.5 0

0.5

1

Fig. 2. Cartoon of the graph of the function Kðu; bÞ of Eq. (7) for five fixed values of b, where b
and bþ are the bifurcation values. The values of b in the subplots are: (a) b ¼ b0 < b
, (b) b ¼ b
, (c) b
< b ¼ b1 < bþ , (d) b ¼ bþ , (e) bþ < b ¼ b2 . Arrows indicate the direction of motion of the system based on the sign of the derivative. Here Iðu; vÞ ¼ v, a ¼ 5, c ¼ 1 and ds ¼ 0:01. Stable critical points are marked with small black circles and unstable points with open circles.

134

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

unstable state. The occurrence of five or more critical points depends on the functional form of I. If the biological system in question is sufficiently complex, such a scenario is conceivable, but our analysis will focus on the one- and two-attractor cases. The stability of a critical point u may be determined by inspecting the relationship between u and b. Implicitly differentiating Eq. (7) with respect to b, the relation     du uoI=ob ¼a ð8Þ db u¼u oK=ou u¼u is obtained. The function I is increasing with b, hence the numerator is positive. The sign of ½du=dbu¼u is therefore equal to the sign of oK=ou in u , i.e. u is unstable whenever ½du=dbu¼u is positive and stable whenever ½du=dbu¼u is negative. 3.3. Bifurcations As outlined earlier, we consider the parameter b to be a measure of the tyrosinase activity, a low (high) value of b corresponds to a low (high) activity of tyrosinase. When b experiences variation, profound changes in the behaviour of Eq. (4a) may take place. Fig. 2(a)–(e) shows a possible scenario: for small values of b there is a unique, stable critical point, when b increases beyond a critical value b
a saddle-node bifurcation occurs as two new critical points emerge, one stable and one unstable. Above b
the system exhibits bistability. When the value of b exceeds a second critical value bþ , the unstable critical point and the stable critical point with the largest value of u merge and disappear and the system returns to monostability. The calculation of the bifurcation values b
and bþ is strongly dependent on the functional form of I. We perform this calculation for the biologically realistic cases when the inactivation of GSH reductase (expressed by the inactivating function I) is (A) assumed to depend on the level of inactivator (v) or (B) on the ratio of inactivator to enzyme concentration (v=u). Moreover, we consider Eq. (4a) in the limit ds ! 0, which corresponds to neglecting the decay of the chemical species glutathione. In this case vðu; bÞ ¼ maxðb
cu; 0Þ and we write, for cases (A) and (B) respectively, I ¼ IA ðvðu; bÞÞ ¼ IA ðmaxðb
cu; 0ÞÞ; I ¼ IB ðvðu; bÞ=uÞ ¼ IB ðmaxðb=u
c; 0ÞÞ; where IA and IB are both increasing with their argument. Setting (A) y ¼ b
cu and (B) y ¼ b=u
c we obtain b ¼ bA ðyÞ ¼ y þ b ¼ bB ðyÞ ¼

c ; 1 þ aIA ðmaxðy; 0ÞÞ

yþc ; 1 þ aIB ðmaxðy; 0ÞÞ

ð9aÞ ð9bÞ

in the equilibrium of Eq. (4a). The interval of bistability is ðb
; bþ Þ on which there are three different equilibrium values of y to each value of b, i.e. there are three critical points of Eq. (4a). The intermediate critical point is unstable. Owing to Eq. (8) and the fact that u and y are inversely related, the interval ðb
; bþ Þ corresponds to the part of the graphs of bA and bB which are decreasing. Hence determining b
and bþ is equivalent to finding the interval ðy
; yþ Þ on which bA

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

135

and bB are decreasing and from the endpoints of this interval, calculating the corresponding values of b. The decreasing part of the graph of bA and bB corresponds to the intermediate, unstable steady state. This state separates the attractor basins of the upper and lower stable steady states. Given the initial state of the system it is possible to assess whether the system is below or above the intermediate state and, from this, determine in which direction the system will move; towards the upper or lower stable state. 3.4. Examples As a first example we consider the inactivating function I ¼ IA ðvÞ ¼ v; which is obtained when the inactivation of GSH reductase is modelled using standard mass action kinetics. We found qualitative correspondence, but no quantitative correspondence with available data on the levels of melanin in murine hairs using this function [17]. In the limit ds ! 0 the inactivation function IA reads IA ðvðu; bÞÞ ¼ maxðb
cu; 0Þ: Letting y ¼ b
cu the function bA ðyÞ (corresponding to the function in Eq. (9a)) reads c : b ¼ bA ðyÞ ¼ y þ 1 þ a maxðy; 0Þ The endpoints of the interval ðy
; yþ Þ on which multiple critical points occur pffiffiffiffiffi is found by locating the local minimum and maximum point of bA . Thus yþ ¼ 0 and y
¼ ð ac
1Þ=a showing that ac > 1 is a necessary requirement for bistability to occur. We mention that this analysis is feasible for positive values of ds . It can be shown that bistability is no longer possible when ds gets large. As a second example we consider the inactivating function I ¼ IB ðv=uÞ ¼

ðv=uÞn vn ¼ ; ðv=uÞn þ d n vn þ ðduÞn

ð10Þ

i.e. the inactivation is assumed to be sigmoidal in v with a threshold du. We found a good quantitative correspondence with available data on the levels of melanin in murine hairs using this function with n ¼ 1:40 [17]. Based on this evidence we consider a model with this functional form to be biologically realistic. In the limit ds ! 0 the function IB reads ½maxð b=u
1; 0Þn ; IB ðvðu; bÞ=uÞ ¼ ½maxð b=u
1; 0Þn þ dn where  b ¼ b=c and d ¼ d=c. Letting y ¼  b=u
1 the function bB ðyÞ (corresponding to the function in Eq. (9b)) reads  b ¼ bB ðyÞ ¼

yþ1 : maxðy; 0Þn 1þa maxðy; 0Þn þ dn

ð11Þ

136

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152 (b)

1.1 b+

1

1

0.8

0.9

0.6

u

β

B

(a)

0.8

0.4

0.7

0.2

b− 0.6

0

10

20

y

30

0

0.6 b

−

0.8

b

1 b 1.2 +

Fig. 3. (a) The function bB ðyÞ in Eq. (11) for the parameter values a ¼ 30, d ¼ 10 and n ¼ 1:40. The points b
and bþ are drawn to illustrate the method of determining the bifurcation values of b. (b) Bifurcation diagram of the system defined by Eqs. (4a) and (5) with inactivation function I of the form (10) and the same parameter values as in (a). The arrows indicate the direction of the quasi-steady state when b increases or decreases.

The endpoints of the interval ðy
; yþ Þ on which multiple critical points occur is found by determining the values of y in which bB attains its local maximum and minimum values. The corresponding b-interval is found by insertion (Fig. 3(a)). The associated bifurcation diagram with the original variable u is given in Fig. 3(b). The intermediate unstable steady state is depicted as the dashed curve in Fig. 3(b). This curve also represents the boundary between the basins of attraction of the coexisting stable states. The equation of the separating curve in terms of b and y is b ¼ bB ðyÞ, where y
< y < yþ . By variable transformation the equation in terms of  b and u can be obtained. 3.5. Bistability in the TMR
We have seen that the TMR0 may exhibit bistability for values of b in the interval ðb
; bþ Þ. Does this mean that the TMR
also has coexisting stable states for b 2 ðb
; bþ Þ? The answer to this question is yes, the main structure of the dynamical system is maintained when the parameter
experiences variation. This is a feature of the specific model being investigated. In general, whenever the small parameter
does not enter the right hand side of the fast subsystem of a singular perturbation problem, the full and reduced systems share the same steady state(s) and therefore the same bifurcations. When
does enter the right hand side of the fast subsystem, this is not necessarily so. In our TMR
it enters in one place, since the second of Eq. (A.6) includes the term a9
½1
ð1 þ aIðxGT ; xDQ ÞÞxGT . Note that the term within the brackets is identical to the right hand side of Eq. (3a). Thus, when we consider the steady state(s) of the TMR
, the term above equals zero. Consequently the steady state(s) of the TMR
and TMR0 are independent of
and

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

137

therefore identical. Furthermore, bistability will emerge at the same bifurcation values in the TMR
as in the TMR0 . It should be emphasised that in general, properties of the TMR
are not necessarily inherited by the corresponding TMR0 . For example, when the smallness parameter
does not disappear when all derivatives are set equal to zero, the steady state(s) of the TMR
will depend on
. Consequently, bifurcations may emerge at different bifurcation values. Hence, in such cases the TMR
can be structurally unstable. 4. Numerical confirmation We complement the development and analysis of the TMR0 by numerical experiments. First, the TMR0 is validated by numerically investigating the convergence properties of the TMR
solution to the solution of the TMR0 . Second, the time evolution of the solution of the TMR0 is considered, given biologically motivated external inputs. 4.1. Comparison and convergence We intend to investigate how rapidly the solution of the TMR
converges towards the solution of the TMR0 . To this end, both models are integrated numerically on the interval ð0; t1 Þ for
k ¼ ð1=2Þk , k ¼ 0; 1; . . . ; 10 for two different values of the parameter ds ; ds ¼ 0:1 and ds ¼ 0:01. The solutions of the two models are depicted in Fig. 4(a) and (b) (details of the numerical experiments are given in the figure legend). In both simulations we observe a marked improvement in the ability of the TMR0 to approximate the solution of the TMR
when
approaches zero, as expected. For each simulation the relative L1 -error ek , i.e. R t1 juðt;
k Þ
xGT ðt;
k Þj dt ek ¼ 0 R t1 ; jxGT ðt;
k Þj dt 0 is computed. The error in this example should be of first order in
, cf. error estimate (B.4). Assuming that ek ¼ C
ak , where C is a constant and a is the order of the error, C and a may be estimated by standard linear regression after taking the logarithm. Using the data for k ¼ 3; . . . ; 9 for ds ¼ 0:1 we obtain the estimates a ¼ 1:0161 and C ¼
0:7966. Similarly, with the data for k ¼ 5; . . . ; 9 when ds ¼ 0:01 the best fit is found for a ¼ 0:7993 and C ¼ 0:2451. We have omitted the data obtained for values of k less than 3 in the linear regression for ds ¼ 0:1 and for values less than 5 in the linear regression for ds ¼ 0:01 due to the slow convergence for large
. The data and the best linear fit are plotted in Fig. 4(c) and (d). In the first simulation, i.e. ds ¼ 0:1, the agreement with error estimate (B.4) is remarkable. In the second simulation, corresponding to ds ¼ 0:01, there is a substantial difference between the numerically obtained value for a ( 0:8) and theory. An inspection of the steady states of the systems with these parameters reveals that for ds ¼ 0:01 the system experiences bistability for values of b in a certain interval. In contrast, for ds ¼ 0:1 the system is monostable regardless the value of b. It therefore appears that the suboptimal rate obtained for ds ¼ 0:01 is related to the bistability of the model for this parameter value, i.e. the convergence of the solution of the TMR
towards the solution of the TMR0 is slightly slower for systems exhibiting bistability than for systems exhibiting monostability. Keeping in mind that the

138

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152 (a) 0.8

1 0.8 xGT

GT

0.6

x

(b)

0.4

0.6 0.4

0.2

0.2

0 0

100

200 Time

300

0 0

400

0

−2

−2

−4 −6 −8 −8

200 Time

300

400

−6

−4 log(ε )

−2

0

(d)

0

log(ek)

log(ek)

(c)

100

−4 −6

−6

−4 log(ε ) k

−2

0

−8 −8

k

Fig. 4. (a,b) Comparison of the numerical solution of xGT found using the TMR
for
¼ 1; 1=2; . . . ; ð1=2Þ10 (thin lines) and the numerical solution of u found using the TMR0 (thicker line) for ds ¼ 0:1 (a) and ds ¼ 0:01 (b). The values of the remaining parameters are: ai , i ¼ 1; . . . ; 11 are all 1, except: a4 ¼ 0:8, a7 ¼ 4:8 and a10 ¼ 1:2. Furthermore, a ¼ 30, c ¼ 1,
0 ¼ 0:1, l ¼ a ¼ 1, hi ¼ 0:5, i ¼ 1; . . . ; 4, the inactivating function I is Iðu; vÞ ¼ v and the time varying levels of aMSH and agouti are, respectively, rm ðt Þ ¼ 5 and ra ðt Þ ¼ 10ð1 þ cosð0:01pt ÞÞ and b1 ðMÞ ¼ M, b2 ðMÞ ¼ b3 ðMÞ ¼ 4M. (c,d) Logarithm of the error ek plotted against the logarithm of
k (o) for ds ¼ 0:1 (c) and ds ¼ 0:01 (d) and the best linear fit to these data.

difference is not dramatic, we conclude that the TMR0 is a trustworthy approximation to the TMR
for small
when the system exhibits monostability as well as when it exhibits bistability. 4.2. Time evolution simulations The model in Eqs. (1) and (2) and the TMR0 (4a) and (5) can be applied to simulate the time evolution of melanin synthesis. Assuming that the levels of external agents a-MSH and agouti protein exhibit temporal variation and imposing this variation on the model should indicate whether the TMR0 is useful for predicting temporal pigment patterns or not. We consider two cases, in each case the level of one of the external agents is varied and one is fixed. In both cases we consider the inactivating function Iðu; vÞ ¼ v corresponding to standard mass action kinetics modelling of the inactivation of GSH reductase. We present the simulation results in terms of the pheomelanin to eumelanin ratio p0 de u ¼c 0 e dp vðu; bÞ þ ds calculated from Eq. (6).

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

139

4.2.1. Variation in agouti protein levels In this simulation the variation in agouti levels is assumed to be Gaussian-like, where the peak may represent a sudden burst in agouti levels occurring midway through the hair growth cycle (see below). We keep the level of a-MSH constant during this simulation. Solving the TMR0 numerically a profile of the ratio of pheomelanin to eumelanin levels (p0 =e0 ) is obtained (Fig. 5(a)). We note that the burst of agouti protein triggers an almost instantaneous jump in the ratio as it switches from a very low level to  2:7 and later drops to the low level. This corresponds to a switching of the system from one basin of attraction to another when the bifurcation parameter exceeds the critical value. The switching may explain the agouti phenomenon in certain mammalian species; hairs of such mammals exhibit a yellow band on a black background, presumably induced by an increase in the expression of agouti half way through the hair growth cycle [12]. 4.2.2. Variation in the a-MSH level We studied the effect of varying levels of a-MSH on melanin production. The level of agouti protein was kept fixed. Solving the TMR0 numerically we obtained a profile of the ratio of pheomelanin to eumelanin levels (Fig. 5(b)). A marked decrease in the p0 =e0 ratio at the time of the onset of a-MSH production was observed, as expected.

(b) 3

2.5

2.5

2

2

p/e

p/e

(a) 3

1.5

1.5

1

1

0.5

0.5

0 2

4

6 Time (days)

8

0 2

4

6 Time (days)

8

Fig. 5. The time dependent ratio of eumelanin to pheomelanin level (solid line) and corresponding (relative) level of 2 agouti protein (a) and a-MSH (b) (dashed line). The levels of external agents are, in (a), rm ðtÞ ¼ 5, ra ðtÞ ¼ 50e
ðt
5:33Þ
ðt
5:33Þ2 and in (b), rm ðtÞ ¼ 50e , ra ðtÞ ¼ 5, where t is the time measured in days after birth. The parameters common to the two simulations are: a ¼ 30, c ¼ 1, ds ¼ 0:01,
0 ¼ 0:01, l ¼ a ¼ 1, b1 ðMÞ ¼ M, b2 ðMÞ ¼ b3 ðMÞ ¼ 4M and de =dp ¼ 0:03.

140

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

5. Discussion 5.1. Bistability behaviour The model of melanogenesis developed by us [17], which has been refined in this paper, is based on the hypothesis that DOPAquinone inactivates the enzyme GSH reductase [16]. This enzyme is crucial for the formation of glutathione (GSH) and thus also for pheomelanogenesis. The statistical treatment in the previous paper [17] and the analytical treatment as well as the numerical simulations of this paper support this basic hypothesis. A good fit to the data of Prota et al. [21] for the genotypic variation of the steady-state level of pheomelanin with tyrosinase activity was obtained [17]. In the present paper we have shown that in an interval of tyrosinase activities, the model may have two stable states. Prior to the development of the TMR0 we were not in a position to analyse the switch between these two states. The analysis of the TMR0 provides an explanation for the eumelanin–pheomelanin switch: at certain bifurcation levels of tyrosinase activity the system switches from one attractor to another. Consequently the melanocyte switches from production of one type of melanin to another over a short range of tyrosinase activity. This is in agreement with our previous findings [17] and the available data [21]. We conclude that our TMR0 really discloses the essential aspects of the melanocyte module. The numerical simulations of Section 4.2 exemplify that the TMR0 may account for and explain the agouti phenomenon seen in some mammals. The abrupt switch may also play a role in the spatial patterning of mammalian coats. The analysis and the numerical simulations indicate that two scenarios may be realised, one where eumelanin dominates, resulting in synthesis of black or dark brown pigment, and one where pheomelanin dominates and the melanocyte produces predominantly yellow or reddish pigment. The spots and stripes in the coat of many mammals are often very sharply delimited. This sharpness may be due to the abrupt nature of melanogenic switching. Thus, even if the agents controlling the production of pigment are relatively smoothly distributed in the epidermal layers, the boundaries between different-coloured regions in the resulting spatial pattern may be very sharp due to the switch-like response of the melanocyte to external input. Given that there are spatial and temporal fluctuations in the external agents and processes affecting the state of a cell and seemingly random variations in processes within the cell, it is reasonable to assume that the state of the cell is stochastic. To cope with this, it is common to infer the probable state of an average cell from samples containing several thousands or millions of cells. Levsky and Singer [22] criticise this methodology as it can only provide information on the Ôaveraged cellÕ, a representation of biological knowledge that does not exist in biology. Does our explanation of melanogenic switching only apply to Ôaveraged melanocytesÕ and not to ÔrealÕ melanocytes? We believe that it does apply to real cells. As discussed above, coat patterns in mammals have very distinct regions of different colour, indicating that the colour determining mechanism is extremely robust towards stochastic variation of external inputs and internal responses. Robustness is a conspicuous feature of the switching mechanism discovered in this paper, supporting the candidacy of this mechanism for the real switching mechanism in melanogenesis. Most models of pattern formation are based on reaction–diffusion mechanisms where patterns emerge due to Turing instabilities e.g. [23,24]. Lacking knowledge on the biological mechanisms

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

141

leading to pattern generation the choice of model is often a choice of convenience rather than a choice based on experimental data. In order to model the formation of coat patterns one could consider the TMR0 to be a module in a larger spatial meta-model of for example the hair follicle system. Using the TMR0 the melanocyte module remains trustworthy and is far more tractable than a module using the full melanogenic model. The TMR0 nicely discloses the output to be expected from varying input signals in accordance with the full model, and in addition it accentuates the main mechanisms determining the relationship between module input and output. However, prior to including the TMR0 in a spatial model the causal relationships between the actors in the TMR0 and in the spatial model have to be determined. 5.2. Normal vs. ‘pathological’ situations In order to obtain the expressions for eumelanin and pheomelanin levels (Eq. (6)) we neglected the degradation of the substances involved in the eumelanin and pheomelanin pathways. In the normal case where the activities of enzymes of these pathways are intermediate or high, degradation may be neglected, but in the ÔpathologicalÕ case where the activities of these enzymes are low, neglecting degradation may cause problems. At low enzyme activities, we can imagine that the conversion proceeds at a lower rate than the degradation, i.e. only a small share of the metabolites is converted to eumelanin or pheomelanin. For example, the activities of the enzymes TRP1 and TRP2 occurring in the eumelanin pathway may vary according to genotype, leading to variation in the level of synthesised eumelanin (see e.g. [21]). This behaviour is not predicted by a model neglecting the degradation of the metabolites. This exposition illustrates that the insights and utility provided by a TMR0 rest on specific assumptions that should not be violated. 5.3. Robustness properties If the TMR0 is to provide valuable information on the behaviour of the biological system, the crucial assumption that the full model correctly describes the system has to be satisfied. Making this assumption, the Levin and Levinson theorem and the error estimate (B.4) ensure that there is a quantitative agreement between the solutions of the TMR0 and the TMR
when
is small. Thus the TMR0 can be considered a trustworthy mathematical representation of the biological system for small
. Unfortunately, there is as yet no clear evidence that the parameter
is small, only indications. In the melanogenesis model, the correspondence is fair between the TMR
and TMR0 solutions for not so small values of
if the external agents vary slowly enough. This suggests that the melanogenic system is remarkably robust with respect to variation in the fast dynamics. This is emphasised by the fact that most of the parameters of the fast subsystem cancel out when the algebraic equations obtained for
¼ 0 are solved. Thus, the making of a TMR0 also represents an efficient sensitivity analysis by analytical means as these parameters appear to be of negligible significance, at least for the qualitative properties of the system. The analysis has also revealed a different kind of robustness in the model. Even though the quantitative aspects of the switching mechanism vary with parameter values, the mechanism itself operates over a wide parameter range (for all
P 0, cf. Section 3.5). In this sense the reduced

142

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

melanogenesis model is a robust module which can readily be incorporated into larger models taking other cell types and spatial extension into account. 5.4. General remarks If
is not small enough to guarantee a quantitative agreement between the TMR0 and the full model, the TMR0 may still provide crucial qualitative insights not so easily seen by inspection of the full model only. If the solution of the TMR0 is close to the solution of the full model in spite of a large
, it is likely that the actual biological system possesses robustness properties. Thus the targeted reduction of a complex model by singular perturbation methods may serve as a means to disclose robustness features that only an extensive numerical investigation and post-simulation analysis of the full model would otherwise be able to provide. However, it should be emphasised that a close correspondence between the TMR0 and the full model does not imply that the system is robust to changes in all fast subprocesses. This is because the algebraic equations representing these processes feed values into the differential equations of the TMR0 . Thus, which type of robustness that is reflected by the correspondence has to be determined in each case, but it is likely that the system is at least robust to perturbations of the parameters determining how fast the subprocesses reach equilibrium after a disturbance. In any case this will be a highly guided search compared to analysis of the full model. Also, if the parameters assumed small really are small, the system will be robust to variations in parameters only affecting these parameters. 5.5. Targeted model reduction TMR
! TMR0 The method applied in this paper is a general procedure applicable to regulatory systems possessing time scale hierarchy. We consider the main steps of a targeted model reduction to be: 1. Setting up a hypothesis on the essential aspects of the model relative to purpose and goal of the modelling. 2. Non-dimensionalising and scaling, exploiting the time scale hierarchy to choose a definite scaling according to the time scale characteristic of the essential processes. This leads to the TMR
. 3. Application of regular and singular perturbation methods to arrive at the TMR0 . 4. Analysis and validation of the TMR0 relative to the real biological system and to the TMR
. Possible benefits are model simplification, improved ways of interpreting and validating the original model, disclosing essential subprocesses and robustness properties, and suggestions as to where experimental efforts should be directed for further validation of the full model. The last aspect includes the fact that the TMR0 is a means to confirm that the correspondence between predictions of the full model and data are really due to the embedded key hypotheses or premises. It is crucial to realise that a specific TMR
version of a large complex model is to a large extent defined according to the given explanatory and utilitarian goals. In this connection the choice of scaling, and especially the time scale, plays an important part. We think that the targeted reduction of complex models will become increasingly important as post-genomic research produces hard biological data on which more and more complex causal

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

143

models can be founded. To get a deep causal understanding of really large systems, access to minimised causal representations of the various subsystems and modules will probably turn out to be quite instrumental.

Appendix A. Model We have recently developed a mathematical model of melanogenesis [17] dealing with binding of a-MSH and agouti on the surface of the melanocyte and with intracellular dynamics. An interaction diagram of intracellular melanogenesis is given in Fig. 1 and abbreviations and mathematical variables are explained in Table 1. The sigmoid Hill function xn ; S ¼ Sðx; hÞ ¼ n x þ hn where the kinetic order n P 1, is used extensively for the modelling of enzyme processes. A.1. Melanogenic pathways within the melanosome This section presents the chemical equations of the subprocesses constituting the process of melanogenesis within the melanosome. The labels of the chemical equations below correspond to the labels given in Fig. 1. More detailed information on the subprocesses are provided in [17]. The two initial steps are common for the eumelanogenic- and the pheomelanogenic pathways and were described by j3

k4

Tyrosine þ Tyrosinase
C1 ! DOPA þ Tyrosinase; k
3

j5

k6

DOPA þ Tyrosinase
C2 ! DQ þ Tyrosinase; k
5

ðC1Þ ðC2Þ

where C1 and C2 are the intermediate enzyme complexes. We considered the major steps of the specific eumelanin pathway to be k7

DQ ! LDC;

ðE1Þ

j8

DQ þ LDC ! DC þ DOPA; j10

ðE2Þ

k11

DC þ TRP2
C3 ! DHICA þ TRP2; k
10

j12

k13

DHICA þ DQ þ TRP1
C4 ! IndQCA þ DOPA þ TRP1; k
12

k14

IndQCA ! Eumelanin; k15

DC ! DHI;

ðE3Þ ðE4Þ ðE5Þ ðE6Þ

144

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152 j16

k17

DHI þ DO
C7 ! IndQ þ DO; k
16

k18

IndQ ! Eumelanin:

ðE7Þ ðE8Þ

Similarly, we considered the major steps of the specific pheomelanin pathway to be j19

DQ þ GSH ! Glu-DOPA; j20

ðP1Þ k21

Glu-DOPA þ c-GTP
C5 ! CDOPA þ c-GTP; k
20

j22

CDOPA þ DQ ! CDQ þ DOPA; k23

CDQ ! Pheomelanin:

ðP2Þ ðP3Þ ðP4Þ

The intermediate GSH is produced from oxidised glutathione (GSSG) by the action of GSHR which we described by the following steps: j24

k25

GSSG þ GSHR
C6 ! GSH þ GSHR: k
24

ðSW1Þ

We assumed that DQ binds itself covalently to GSHR and thereby reduces GSHR activity, i.e. j26

DQ þ GSHR
GSHRB ; k
26

ðSW2Þ

where GSHRB denotes inactivated GSHR. A.2. Dimensional model The model of the dynamics on the surface of the melanocyte was described by _ ¼ j1 Rm ðtÞð1
M
AÞ
k
1 M; M A_ ¼ j2 Ra ðtÞð1
M
AÞ
k
2 A:

ðA:1Þ

We modelled the first steps of intramelanosomal melanogenesis by T_ ¼ s
cT ; D_ p ¼ hj3 ST T
/j5 Dp T þ j8 Dq Lc þ SðDq Da ; Kqa Þk13 ETRP 1 þ j22 Dq Cp ; D_ q ¼ /j5 Dp T
k7 Dq
j8 Dq Lc
SðDq Da ; Kqa Þk13 ETRP 1
j19 Dq Gs
j22 Dq Cp þ ½k
26 ðEGSHR
GT Þ
k~26eI ðGT ; Dq ÞGT ;

ðA:2Þ

where h ¼ k4 =ðk
3 þ k4 Þ and / ¼ k6 =ðk
5 þ k6 Þ and the function I~ describes the binding of DQ to GSH reductase. We assumed that the tyrosine concentration ST is constant and that the parameters reflecting tyrosinase activity and synthesis, j3 , j5 and s, are functions of the fractional occupancy of MC1-R by a-MSH (M).

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

145

The eumelanin pathway was modelled by L_ c ¼ k7 Dq
j8 Dq Lc ; D_ c ¼ j8 Dq Lc
SðDc ; Kc Þk11 ETRP 2
k15 Dc ; D_ a ¼ SðDc ; Kc Þk11 ETRP 2
SðDq Da ; Kqa Þk13 ETRP 1 ; I_a ¼ SðDq Da ; Kqa Þk13 ETRP 1
k14 Ia ; D_ h ¼ k15 Dc
SðDh ; Kh Þk17 EDO ; I_ ¼ SðDh ; Kh Þk17 EDO
k18 I;

ðA:3Þ

E_ ¼ k14 Ia þ k18 I
cE E: Finally, the pheomelanin pathway was described by G_ T ¼ k
26 ðEGSHR
GT Þ
k~26 IðGT ; Dq ÞGT ; G_ s ¼ k~25 GT
j19 Dq Gs
cs Gs ; G_ p ¼ j19 Dq Gs
SðGp ; Kp Þk21 Ec-GTP ; C_ p ¼ SðGp ; Kp Þk21 Ec-GTP
j22 Dq Cp ; C_ q ¼ j22 Dq Cp
k23 Cq ; P_ ¼ k23 Cq
cP P :

ðA:4Þ

A.3. Non-dimensional model For the non-dimensionalisation of the model we define the slow dimensionless time t ¼ k
26 t. The dimensionless variables are defined in Table 1 and the parameters in Table 2. The subsequent non-dimensional model is given by Eqs. (1), (2), (3a) and (3b). Here we define the variables y and the right hand side g occurring in Eq. (3b). We set y ¼ ðxc ; xe ; xp Þ. The vectors xc , xe and xp represent the variables associated with the common pathway and the eumelanin and pheomelanin pathways, respectively. The variables xe and xp are defined in Eqs. (A.7) and (A.9) below, respectively, and xc ¼ ðxDP ; xDQ Þ:

ðA:5Þ

Similarly, we set g ¼ ðgc ; ge ; gp Þ, where the vector functions gc , ge and gp incorporate the interactions between the metabolites of the common pathway and the eumelanin and pheomelanin pathways, respectively, with other metabolites. These vector functions are defined in Eqs. (A.6), (A.8) and (A.10) below, respectively. The components gcð1Þ and gcð2Þ of gc are gcð1Þ ¼ b2 ðMÞxT
b3 ðMÞxT xDP þ xDQ xLDC þ a1 SðxDQ xDA ; h21 Þ þ a2 xDQ xCP ; gcð2Þ ¼ b3 ðMÞxT xDP
xDQ
xDQ xLDC
a1 SðxDQ xDA ; h21 Þ
xDQ xGS
a2 xDQ xCP þ a9
½1
xGT
aIðxGT ; xDQ ÞxGT ;

ðA:6Þ

146

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

where IðxGT ; xDQ Þ ¼ eI ðGT ; Dq Þ. The vector xe is xe ¼ ðxLDC ; xDC ; xDA ; xIA ; xDH ; xI ; eÞ

ðA:7Þ

and the components of the function ge read geð1Þ ¼ a3 ð1
xLDC ÞxDQ ; geð2Þ ¼ xDQ xLDC
a4 SðxDC ; h2 Þ
a5 xDC ; geð3Þ ¼ a4 SðxDC ; h2 Þ
a1 SðxDQ xDA ; h21 Þ; geð4Þ ¼ a1 SðxDQ xDA ; h21 Þ
a6 xIA ; geð5Þ geð6Þ geð7Þ

ðA:8Þ

¼ a5 xDC
a7 SðxDH ; h3 Þ; ¼ a7 SðxDH ; h3 Þ
a8 xI ; ¼ a6 xIA þ a8 xI
de e:

The vector xp is xp ¼ ðxGS ; xGP ; xCP ; xCQ ; pÞ

ðA:9Þ

and the components of the function gp read gpð1Þ ¼ cxGT
xDQ xGS
ds xGS ; gpð2Þ ¼ xDQ xGS
a10 SðxGP ; h4 Þ; gpð3Þ ¼ a10 SðxGP ; h4 Þ
a2 xDQ xCP ;

ðA:10Þ

gpð4Þ ¼ a2 xDQ xCP
a11 xCQ ; gpð5Þ ¼ a11 xCQ
dp p:

Appendix B. Mathematical background B.1. Time scales in models of biological systems Frequently, models of biological systems comprise subprocesses with typical time scales ranging over many orders of magnitude. In the following we outline briefly how this feature may be employed to arrive at a new representation of the model more amenable to analysis than the original. We emphasise that this appendix is intended to be a means of tuning the reader unfamiliar with singular perturbation methods into the right frequency. Readers who are familiar with singular perturbation theory may well skim through or entirely skip the appendix. The basis for the analysis is a hypothesis that certain of the processes in the model are essential for those aspects of the real system that the model is intended to grasp, while other subprocesses are of less interest. Typically realistic models for biological systems comprise large numbers of parameters, the majority of which are either completely unknown or subject to rough estimates of magnitude. For simplicity we let the parameters be sorted according to their assumed magnitude in just two sets, corresponding to two typical time scales, i.e. subprocesses in the model occur roughly either fast or at rates comparable to the essential processes of the model.

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

147

Consider a real system for which there exists a non-dimensionalised model dX ¼ Uðt; X ; P ; QÞ; ðB:1Þ dt where X 2 RN , and P and Q are the two sets of parameters. The parameters in these sets typically represent relative rates stemming from the biological system on which model (B.1) is based. Moreover, we assume that the parameters in P and Q can be associated with time scales of the system. Parameters that cannot be associated with time scales (e.g. Hill coefficients) are not of interest and are therefore omitted in Eq. (B.1). Let Pi and Qj be any elements of P and Q and assume that Pi  Qj for all i; j, while all elements of each of the two sets are roughly equal in size. The biological/biochemical interpretation of this is that the parameters in P and Q represent fast and slow processes, respectively. To exploit this fact we select two parameters p and q among the elements P and Q, respectively, and introduce a small parameter
¼ q=p. The idea is that q is a relative rate representing the time scale of what we believe are the key regulatory processes of the system. Thus, we also choose a time scale accordingly and define a scaled time t by t ¼ qt . To simplify the presentation assume that there is just a single parameter p in P . Then p ¼ q=
and dX ¼ Uðt; X ; q=
; QÞ: dt We keep q fixed and consider what happens when
approaches zero. Normally the components of U will either stay non-zero, but bounded, or tend to infinity. Thus, we group the components of U in two functions F and G according to this behaviour; F represents the components of U which are bounded when
! 0 and G represents the components which tend to infinity. Assume for simplicity that for 0 <
 1 kF ðt; X ; q=
; QÞk < MF ; 1 kGðt; X ; q=
; QÞk MG ;
where MF and MG are positive constants and k  k defines the usual Euclidean norm. Let x and y be the variables whose rate functions are F and G, respectively. We define the functions f and g by f ðt ; x; y;
; QÞ ¼ qF ðt; X ; q=
; QÞ; gðt ; x; y;
; QÞ ¼
qGðt; X ; q=
; QÞ: These functions are non-zero and bounded when
approaches zero. Including initial conditions, this leads to the equivalent model representation dx ¼ f ðt ; x; y;
; QÞ; xðt0 Þ ¼ a; dt ðB:2Þ dy
 ¼ gðt ; x; y;
; QÞ; yðt0 Þ ¼ b: dt In this representation the model is ready for regular and singular perturbation analysis, for which a large body of theory exists (see Appendix B.2 below). In singular perturbation literature, system (B.2) is commonly referred to as the full problem. The first equation of (B.2) represents the normal subsystem and the second the rapidly moving subsystem. If
is small, one would intuitively expect an approximate solution to be given by the

148

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

reduced problem (some authors also refer to it as the degenerate problem), the algebro-differential initial-value problem dx0 ¼ f ðt ; x0 ; Wðt ; x0 Þ; 0; QÞ; dt y0 ¼ Wðt ; x0 Þ;

xðt0 Þ ¼ a;

ðB:3Þ

where the singular state Wðt ; x0 Þ is obtained by solving gðt ; x0 ; y0 ; 0; QÞ ¼ 0 and x0 and y0 are vectors of the same dimensions as x and y. B.2. Convergence theorems Based on assumptions concerning existence and uniqueness of solutions, smoothness of rate functions and stability of the fast subsystem, Friedrich and Wasow [25] proved convergence results for system (B.2) in relation to system (B.3) when
! 0 valid when x is a vector of arbitrary dimension, but y is of dimension one. Their stability assumption was of the simple form og=oy 6
m < 0 when ðx; yÞ is restricted to a certain set and m > 0. Levin and Levinson [26] studied the generalisation that y is a vector of arbitrary dimension. The assumptions on which their convergence result (see below) was based, were the following. Assumption 1. The reduced system possesses a solution ðx; yÞ ¼ ðx0 ðtÞ; y0 ðtÞÞ which is continuously differentiable on the closed and finite interval t0 6 t 6 T0 . Assumption 2. The functions f and g are continuously differentiable with respect to ðt; x; y;
Þ in a region R defined by R ¼ fðt; x; y;
Þ : kx
x0 ðtÞk < d0 ; ky
y0 ðtÞk < d0 ; t0 6 t 6 T0 ; 0 6
< d0 g; where d0 > 0. Assumption 3. For all eigenvalues k of the Jacobian corresponding to the fast subsystem, the inequality Re k 6
l < 0 holds, where l > 0. The Levin and Levinson theorem [26] is formulated as follows. Theorem (Theorem 1:2 in [26]). Assume the initial value problems given by systems (B.2) and (B.3) satisfy Assumptions 1–3. If
is sufficiently small, then system (B.2) has a unique solution ðxðt;
Þ; yðt;
ÞÞ existing over the whole interval t0 6 t 6 T0 and kxðt;
Þ
x0 ðtÞk þ kyðt;
Þ
y0 ðtÞk ! 0 uniformly over t1 6 t 6 T0 as
! 0, for any t1 > t0 . The practical importance of the theorem is that for small
and t > t0 the solution of the reduced problem (B.3) is a good approximation to the solution of the full problem (B.2). If this is the case

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

149

the solution of the reduced problem is a uniformly valid approximation. The Levin and Levinson Theorem is crucial for the development and validity of our TMR0 . The accuracy of the solution of the reduced initial value problem is described by the following estimate [27, Section 6.3], where ðx; yÞ and ðx0 ; y0 Þ are the solutions of problems (B.2) and (B.3), respectively: kxðt;
Þ
x0 ðtÞk; kyðt;
Þ
y0 ðtÞk 6 C
:

ðB:4Þ

Here C > 0 is a constant. The estimate holds for t > 0 and
6
0 where
0 is a suitable, fixed, positive constant. The convergence theorem most commonly cited for problems of the type (B.2) and (B.3) is probably that of Tikhonov [28] (see e.g. [29, Section 39], for an English translation). With a set of stricter assumptions than Levin and Levinson and Tikhonov, Hoppensteadt [30,31] showed the important generalisation that the convergence is uniform for infinitely wide t-intervals. Appendix C. Proof of the main result C.1. Application of Levin and LevinsonÕs theorem We intend to show that the solution of the TMR
(Eqs. (3a) and (3b)) converges uniformly towards the solution of the TMR0 (Eqs. (4a) and (5)) as
! 0, i.e. we intend to show that the TMR0 is a uniformly valid approximation to the TMR
. We start with a simple observation. Remark. By inspection the right hand sides of Eqs. (3a) and (3b) are continuous and differentiable (see Eqs. (A.6), (A.8), and (A.10)), asserting that Assumption 2 of the Levin and Levinson theorem is fulfilled. Assumption 3 is the crucial stability assumption of the Levin and Levinson theorem. Let J be the Jacobian of system (3b), where the matrix elements are the derivatives of g with respect to the components of y evaluated in y ¼ y0 . Unfortunately, we have not been able to prove that J satisfies Assumption 3. However, Monte Carlo simulations show that the real part of all eigenvalues of J are bounded above by a negative value (see Appendix C.2), leading to the following. Conjecture. The real part of all eigenvalues of J are negative, i.e. for all eigenvalues k the inequality Rek 6
l < 0 holds, where l > 0. To assert that Assumption 1 of the Levin and Levinson theorem is fulfilled a result on the regularity of the solution of the reduced problem is necessary. Lemma. The solution of the TMR0 Eqs. (4a) and (4b) is continuous and differentiable for all t . Proof. Observe that the solution of Eq. (4b), y0 , is a continuous and differentiable function of u. In particular v (the second component of y0 ) is a continuous and differentiable function of u. Thus, Iðu; vðu; bÞÞ is also a continuous and differentiable function of u. In order to show that the solution

150

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

of the TMR0 is continuous and differentiable with respect to t we only need to consider the solution of Eq. (4a). The corresponding integral equation is Z t  ½1
½1 þ aIðuðsÞ; vðuðsÞ; bÞÞuðsÞ ds: uðt Þ ¼ uðt0 Þ þ t0

Thus, for arbitrary t1 and t2 ,   Z t2    juðt2 Þ
uðt1 Þj ¼  ½1
½1 þ aIðuðsÞ; vðuðsÞ; bÞÞuðsÞ ds 6 ð1 þ aIð0; bÞÞjt2
t1 j; t1

which shows that u is continuous for all t . Furthermore,   du  du 
  ¼
ðuðt2 Þ
uðt1 ÞÞ
a½bI ðuðt2 ÞÞ
bI ðuðt1 ÞÞ; dt  dt t¼t2

t¼t1

where bI ðuÞ ¼ uIðu; vðu; bÞÞ is a continuous function of u. By the continuity of bI and u, the right hand side approaches zero when t2 approaches t1 . Hence u is differentiable. This concludes the proof. h The Remark, the Lemma and the Levin and Levinson theorem now give the main result in Section 2.3. C.2. On the eigenvalues of the Jacobian Upon elimination of the rows and columns associated with real reduced Jacobian of the function g of Eq. (3b) reads 8
A 1 þ B þ xGS 0 D 0 > > > A
ð2 þ B þ 2x Þ  > 0
D
x > GS DQ > > > 0 1
ðC þ a5 Þ 0 0 < e 0
B C
D 0 J ¼ >   > 0 0
ðxDQ þ ds Þ 0
xGS > > >  > 0 0 xDQ 0 x > GS > :  0 0 0 0
xGS

and negative eigenvalues the 0 0 0 0 0
E E

9 a2 xDQ > > >
a2 xDQ > > > > 0 > = 0 ; > 0 > > > > 0 > > >  ;
a2 xDQ

where the parameters denoted by uppercase letters are A ¼ b3 ðMÞxT ;

B ¼ na1 S1 ð1
S1 Þ=xDQ ;

D ¼ na1 S1 ð1
S1 Þ=xDA ;

C ¼ na4 S2 ð1
S2 Þ=xDC ;

E ¼ na10 S3 ð1
S3 Þ=xGP

and n

S1 ¼

ðxDQ xDA Þ

ðxDQ xDA Þn þ ðh21 Þ

n;

S2 ¼

ðxDC Þn ; n ðxDC Þ þ hn2

S3 ¼

ðxGP Þn : n ðxGP Þ þ hn4

Our Monte Carlo simulations proceed as follows. The quantities xGS , b2 ðMÞxT , a2 , a5 , A, B, C, D and E appearing in the Jacobian are randomly chosen from the interval ½K0 ; K1  and ds is picked at random from ½d0 ; d1 . (The numerical values of K0 , K1 , d0 and d1 are given in the legend of Fig. 6). J with the chosen entries are computed. Moreover, xDQ ¼ b2 ðMÞxT =ð1 þ xGS Þ. The eigenvalues of e

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

x 10

−3

Qk

−1

151

−2

0

500

1000 k

1500

2000

Fig. 6. The elements Qk as a function of the index k. Here, K0 ¼ 0:1, K1 ¼ 10, d0 ¼ 0:001 and d1 ¼ 0:1. Furthermore, L ¼ 1000 and N ¼ 2000.

The eigenvalue kk;j with the largest real value is determined and the real part of this eigenvalue L pk;j ¼ Rekk;j is stored. This procedure is repeated for j ¼ 1; . . . ; L to produce the sequence fpk;j gj¼1 from which we determine Pk ¼ max1 6 j 6 L pk;j , the largest real value of any of the eigenvalues computed in simulation k. This simulation is performed for k ¼ 1; . . . ; N resulting in the sequence fPk gNk¼1 . In order to investigate whether there is an upper limit on the real part of the eigenvalue we define the increasing sequence fQk gNk¼1 by 1. If k ¼ 1: Q1 ¼ P1 , 2. If k > 1: Qk ¼ maxðPk ; Qk
1 Þ. If fQk g is increasing and converging to an upper bound Q < 0 when k ! N , this would support the conjecture. The plot of Qk as a function of the index k (Fig. 6) clearly indicates that the eigenvalues all have negative real part and that the limit Q ¼ limk!1 Qk exists and Q < 0. To keep the number of parameters low we chose parameters A, B, C, D and E at random despite the fact that these parameters depend on xGS and xDQ . Thus the result of the simulations holds for a larger class of matrices than necessary. This could mean that only a very small portion of the simulations describe the situation of interest to us, and thereby weaken the conclusion. However, having repeated the simulations 2 · 106 times, we consider it unlikely that our conclusion is wrong. References [1] J.D. Murray, Mathematical Biology, Springer, Heidelberg, 1983. [2] J.G. Reich, E.E. SelÕkov, Mathematical analysis of metabolic networks, FEES Lett. 40 (Suppl.) (1974) S119.

152

L. Øyehaug et al. / Mathematical Biosciences 185 (2003) 123–152

[3] J.G. Reich, E.E. SelÕkov, Time hierarchy equilibrium and non-equilibrium in metabolic systems, Biosystems 7 (1974) 39. [4] R. Heinrich, I. Sonntag, Dynamics of non-linear biochemical systems and the evolutionary significance of time hierarchy, Biosystems 15 (1982) 301. [5] L.A. Segel, On the validity of the steady state assumption of enzyme kinetics, Bull. Math. Biol. 50 (1988) 579. [6] F. Battelli, C. Lazzari, Singular perturbation theory for open enzyme reaction networks, IMA J. Math. Appl. Med. Biol. 3 (1986) 41. iska, Analysis of metabolic systems with complex slow and fast dynamics, Bull. Math. Biol. 51 [7] I. Dvorak, J. S (1989) 255. [8] J.C. Liao, E.N. Lightfoot, Extending the quasi-steady state concept to analysis of metabolic networks, J. Theor. Biol. 126 (1987) 253. [9] J. Delgado, J.C. Liao, Control of metabolic pathways by time-scale separation, Biosystems 36 (1995) 55. [10] L.W. Maki, J. Keizer, Mathematical analysis of a proposed mechanism for oscillatory insulin secretion in perifused HIT-15 cells, Bull. Math. Biol. 57 (1995) 569. [11] Y. Tang, J.L. Stephenson, H.G. Othmer, Simplification and analysis of models of calcium dynamics based on IP3sensitive calcium channel kinetics, Biophys. J. 70 (1996) 246. [12] S.J. Bultman, E.J. Michaud, R.P. Woychik, Molecular characterization of the mouse agouti locus, Cell 71 (1992) 1195. [13] D.S. Lu, D. Willard, I.R. Patel, S. Kadwell, L. Overton, T. Kost, M. Luther, W.B. Chen, R.P. Woychik, W.O. Wilkison, R.D. Cone, Agouti protein is an antagonist of the melanocyte-stimulating-hormone receptor, Nature 371 (1994) 799. [14] K.G. Mountjoy, L.S. Robbins, M.T. Mortrud, R.D. Cone, The cloning of a family of genes that encode the melanocortin receptors, Science 257 (1992) 1248. [15] L.S. Robbins, J.H. Nadeau, K.R. Johnson, M.A. Kelly, L. Rosel-lirehfuss, E. Baack, K.G. Mountjoy, R.D. Cone, Pigmentation phenotypes of variant extension locus alleles result from point mutations that alter MSH receptor function, Cell 72 (1993) 827. [16] S. Ito, Biochemistry and Physiology of Melanin, in: N. Levine (Ed.), Pigmentation and Pigmentary Disorders, CRC, Boca Raton, FL, 1993, p. 34. [17] L. Øyehaug, E. Plahte, D.I. V
age, S.W. Omholt, The regulatory basis of melanogenic switching, J. Theor. Biol. 215 (2002) 449. [18] R. Thomas, Logical vs continuous description of systems comprising feedback loops: the relation between time delays and parameters, Stud. Phys. Theor. Chem. 28 (1983) 307. [19] E.H. Snoussi, R. Thomas, Logical identification of all steady states: The concept of feedback loop characteristic states, Bull. Math. Biol. 55 (1993) 973. [20] E. Plahte, T. Mestl, S.W. Omholt, Feedback loops, stability and multistationarity in dynamical systems, J. Biol. Syst. 3 (1995) 409. [21] G. Prota, M.L. Lamoreux, J. Muller, T. Kobayashi, A. Napolitano, M.R. Vincensi, C. Sakai, V.J. Hearing, Comparative-analysis of melanins and melanosomes produced by various coat color mutants, Pigment Cell Res. 8 (1995) 153. [22] J.M. Levsky, R.H. Singer, Gene expression and the myth of the average cell, Trends Cell Biology 13 (2002) 4. [23] J.D. Murray, On pattern formation mechanisms for lepidopteran wing patterns and mammalian coat markings, Philos. Trans. R. Soc. Lond. B. Biol. Sci. 295 (1981) 473. [24] J.D. Murray, How the leopard got its spots, Sci. Am. 258 (1988) 80. [25] K.O. Friedrichs, W.R. Wasow, Singular perturbation of non-linear oscillations, Duke Math. J. 13 (1946) 367. [26] J.J. Levin, N. Levinson, Singular perturbations of non-linear systems of differential equations and an associated boundary layer equation, J. Rational Mech. Anal. 3 (1954) 247. [27] D.R. Smith, Singular-Perturbation Theory: An Introduction with Applications, Cambridge University, Cambridge, 1985. [28] A.N. Tikhonov, Systems of differential equations containing a small parameter multiplying the derivative, Mat. Sb. N. S. 73 (5) (1952) 575 (In Russian). [29] W. Wasow, Asymptotic Expansions for Ordinary Differential Equations, Dover, New York, 1987. [30] F. Hoppensteadt, Singular perturbations on the infinite interval, Trans. Am. Math. Soc. 123 (1966) 521. [31] F. Hoppensteadt, Asymptotic stability in singular perturbation problems, J. Diff. Eqs. 4 (1968) 350.