## Abstract

Mathematical models are increasingly being used to understand complex biochemical systems, to analyse experimental data and make predictions about unobserved quantities. However, we rarely know how robust our conclusions are with respect to the choice and uncertainties of the model. Using algebraic techniques, we study systematically the effects of intermediate, or transient, species in biochemical systems and provide a simple, yet rigorous mathematical classification of all models obtained from a *core model* by including intermediates. Main examples include enzymatic and post-translational modification systems, where intermediates often are considered insignificant and neglected in a model, or they are not included because we are unaware of their existence. All possible models obtained from the core model are classified into a finite number of classes. Each class is defined by a mathematically simple *canonical model* that characterizes crucial dynamical properties, such as mono- and multistationarity and stability of steady states, of all models in the class. We show that if the core model does not have conservation laws, then the introduction of intermediates does not change the steady-state concentrations of the species in the core model, after suitable matching of parameters. Importantly, our results provide guidelines to the modeller in choosing between models and in distinguishing their properties. Further, our work provides a formal way of comparing models that share a common skeleton.

## 1. Introduction

Systems biology aims to understand complex systems and to build mathematical models that are useful for inference and prediction. However, model building is rarely straightforward, and we typically seek a compromise between the simple and the accurate, shaped by our current knowledge of the system. Two models of the same system, potentially differing in the number of species and the form of reactions, might have different qualitative properties and the conclusions we draw from analysing the models might be strongly model-dependent. The predictive value and biological validity of the conclusions might thus be questioned. It is therefore important to understand the role and consequences of model choice and model uncertainty in modelling biochemical systems.

Transient, or intermediate, species in biochemical reaction pathways are often ignored in models or grouped into a single or few components, either for reasons of simplicity or conceptual clarification, or because of lack of knowledge. For example, models of the multiple phosphorylation systems vary considerably in the details of intermediates [1,2], and intermediates are often ignored in models of phosphorelays and two-component systems [3,4]. Typically, intermediate species are protein complexes such as a kinase–substrate protein complex. It has been shown that sequestration of intermediates can cause ultrasensitive behaviour in some systems [5,6]. Therefore, the inclusion/exclusion of intermediates is a matter of considerable concern.

As an example, consider the transfer of a modifier molecule, such as a phosphate group in a two-component system, from one molecule to another: , where *A*, *B* are unmodified forms (without the modifier group), *A**, *B** are modified forms (with the modifier), indicate reversible reactions, and *…* are potential transient reaction steps. Two-component systems are ubiquitous in nature and vary considerably in architecture and mechanistic details across species and functionality [7]. Whether or not the specifics are known beforehand, it is customary to use a reduced scheme such as [3,4].

We use chemical reaction network theory (CRNT) to model a system of biochemical reactions and assume that the reaction rates follow mass-action kinetics. The polynomial form of the reaction rates has made it possible to apply algebraic techniques to learn about qualitative properties of models, without resorting to numerical approaches [8–14]. Building on previous works [15–17], we propose a mathematical framework to compare different models and to study the dynamical properties of models that differ in how intermediates are included. The most fundamental and crucial dynamical features are the number and stability of steady states. We assume that the kinetic parameters are unknown and study the capacity of each model to exhibit different steady-state features.

The paper is organized in the following way. We first introduce the concepts of a core model and an extension model. An extension model is constructed from the core model by including intermediates. Next, we discuss how the steady-state equations of different models are related and illustrate the findings with an example. We proceed to discuss the number of steady states of core and extension models. After that, we introduce the steady state classes, a key concept of this paper. Extension models in the same steady-state classes have the same properties at steady state (provided the parameter sets of the two models can be matched, in some sense). Using these ideas, we build a decision tree to guide the modeller in choosing a model and in understanding the consequences of choosing a particular model. Finally, we illustrate our approach with an example based on two-component systems. All proofs and mathematical details are in the electronic supplementary material.

## 2. The core model and its extensions

We use the notation and formalism of CRNT [18,19]. A *reaction network* is defined as a set of species, denoted by capital letters (e.g. *A*, *B*, *C*), a set of complexes and a set of reactions between complexes. Each complex is a combination of species, for example *y*_{1} = *A* + *B* or *y*_{2} = 2*C* (not to be confused with a protein complex). A potential reaction could be *A* + *B →* 2*C*, or also written simply *y*_{1} → *y*_{2}. A reaction is not necessarily reversible, that is, we can have *A* + *B →* 2*C* without having the reverse reaction 2*C → A* + *B*. Whenever a reaction is reversible, we model it as two separate irreversible reactions. We assume that each reaction occurs according to mass-action kinetics, that is, at a rate proportional to the product of the species concentrations in the reactant or source complex [20]. For example, the reaction *A* + *B →* 2*C* occurs at a rate *k*[*A*][*B*], where [*A*], [*B*] are the concentrations of the species *A*, *B* and *k* is a reaction-specific positive constant. Reaction networks are often drawn graphically as shown in figure 1*a–e*. Figure 1*c–e* is schematic of reaction networks: only the structure of the network is shown, and neither the species nor the rate constants are indicated.

Figure 1*a* corresponds to a simple enzymatic mechanism where *E* is an enzyme and *S _{j}* is a substrate with

*j*= 0,1,2 phosphorylated sites. The substrate

*S*

_{0}can be doubly phosphorylated sequentially via

*S*

_{1}or directly (processively). In figure 1

*b*, a transient product

*Y*formed by

*S*

_{0}and

*E*, or by

*S*

_{1}and

*E*(these are often denoted by

*S*

_{0}

*· E*and

*S*

_{1}

*· E*) is shown. In the particular case, we do not distinguish between the two transient products (which might be unrealistic, but it serves an illustrative purpose).

An *intermediate* is defined as a species in a reaction network that is created and dissociated in isolation, that is, it is produced in at least one reaction, consumed in at least one reaction and it cannot be part of any other complex (e.g. *Y* in figure 1*b–d*). A *core model* is the minimal reaction mechanism to be modelled. Each reaction *y _{i} → y_{j}* in the core model consists of two

*core complexes y*,

_{i}*y*. The species contributing to the core complexes are referred to as

_{j}*core species*. An

*extension model*is any reaction network such that

(i) the set of complexes consists of core complexes and some intermediates that are not part of the core model;

(ii) reactions are between two core complexes, two intermediates or between an intermediate and a core complex; and

(iii) the core model is obtained from the extension model by

*collapsing*all reaction paths*y*_{i}→*Y*_{1}→*…*→*Y*→_{k}*y*, where_{j}*Y*are intermediates, into a single reaction_{i}*y*→_{i}*y*._{j}

Some examples are given in figure 1. Figure 1*b* is an extension model of figure 1*a* and figure 1*c* and *d* are extension models of figure 1*e*. Figure 1*a* is a concretization of figure 1*e*. Note that the directionality of the reaction arrows needs to be preserved. For instance, in figure 1*e*, an extension of the reaction *y*_{1} → *y*_{2} cannot be , because it would imply that *y*_{2} → *y*_{1} also is in the core model. By adding arbitrarily many intermediates (e.g. ), we can create arbitrarily many extension models with the same core.

Under mass-action kinetics, the dynamics of figure 1*b* is described by a polynomial system of ordinary differential equations (ODEs):
2.1where *k*_{*} are rate constants, [*X*] denotes the concentration of species *X* and is the instantaneous change in [*X*]. In addition, there are two *conservation laws*,
2.2that is, quantities that are conserved over time and determined by the initial concentrations. Conservation laws confine the dynamics to an invariant space given by and (referred to as *conserved amounts*), and the dynamical analysis must be restricted to this space. The invariant spaces are called *stoichiometric classes* in the CRNT literature. If we consider a maximal set of independent conservation laws, then the species that appear in the conservation laws are independent of the chosen set.

The core model in figure 1*a* has two conservation laws,
2.3The two sets of conservation laws, (2.2) and (2.3), differ by a linear combination of intermediate concentrations (here a single term). This similarity between (2.2) and (2.3) holds generally.

### Theorem 2.1.

*The conservation laws in the core model are in one-to-one correspondence with the conservation laws in any extension model. The correspondence is obtained by adding a suitable linear combination of the [Y]'s to each conservation law of the core model.*

Theorem 2.1 does not depend on the assumption of mass-action kinetics but relies on the structure of the network only, that is, on the set of reactions of the network.

## 3. Steady-state equations

We next state theorems 3.1 and 3.2 that allow us to relate the dynamics near steady states of the core and extension models to each other.

At steady state for all species *X*. Under the assumption of mass-action kinetics, this condition translates into a system of polynomial equations in the species concentrations. A way to solve the equations is to express one variable in terms of other variables. This expression must then be satisfied by any solution to the system. We let [*y*] denote the product of the species concentrations in complex *y*, for example, [2*S*] = [*S*]^{2} and [*S*_{0} + *E*] = [*S*_{0}][*E*]. Different extension models contain different intermediates, resulting in different steady-state equations. Because the intermediates always appear as linear terms in the steady-state equations of an extension model (see (2.1)), they can be eliminated from the equations and written in terms of the concentrations of the core species.

### Theorem 3.1.

*Using the equations* *for all intermediates in the extension model, the steady-state concentrations of the intermediates Y are given as linear sums* *of products of the core species concentration* [17,21]*. The constant μ_{Y,y} is either zero or positive and depends only on the rate constants of the extension model. [y] appears in the expression, that is, μ_{Y,y} ≠ 0, if and only if there is a reaction path y→ … →Y involving exclusively intermediates*.

As a consequence of the theorem, once the steady-state concentrations of the core species are known, the steady-state concentrations of the intermediates are also known. Because *μ*_{Y,y} ≥ 0 and at least one of the constants is non-zero (all intermediates are produced), positive steady-state concentrations of the core species lead to positive concentrations of the intermediates.

Theorem 3.1 makes explicit use of mass-action kinetics. It remains true for non-mass action kinetics in the sense that an explicit expression for [*Y*] can be found if all reactions have mass-action reaction rates, whereas all other reactions can have arbitrary reaction rates. In that case, however, the form of the expression might not be polynomial nor lead to positive concentrations.

The manipulations leading to the expression from are purely algebraic and do not require any assumptions about the conserved amounts. In example (2.1), the equation gives
3.1where are reciprocal Michaelis–Menten constants [20]. If (3.1) is substituted into (2.1), we obtain a new ODE system:
3.2which is a mass-action system for the core model in figure 1*a* with , and (as ). We say that the rate constants are *realized* by and that *k*_{*} and are a pair of *matching* rate constants. In the particular case, are realized by choosing , and any , such that . Choosing *k*_{2} fixes the values of *m*_{1}, *m*_{3} in (3.1). However, for some (unrealistic) extension models, not all choices of rate constants of the core model are realizable (see electronic supplementary material).

The relation between the ODEs in figure 1*a,b* holds generally for any pair of core and extension models.

### Theorem 3.2.

*After substituting the expressions* *into the ODEs of the extension model, we obtain a mass-action system for the core model*.

The quasi-steady-state approximation (QSSA) proceeds similarly [22]. An equation of the form is used to find an expression for [*Y*] in terms of [*y*] under the additional assumptions that certain species are in high or low concentration. This expression is subsequently substituted into the remaining ODEs to reduce the system. Theorems 3.1 and 3.2 show that this always can be performed, irrespective of any biological justification of the procedure.

As a consequence of the theorems, the steady states of an extension model can be found in this way: we first solve the equations for [*Y*] in terms of [*y*] (theorem 3.1) and then insert the expressions for [*Y*] into the remaining steady-state equations (theorem 3.2). The steady states of the extended model are now found by solving the steady-state equations for the core model to obtain the concentrations of the core species. This corresponds to solving (3.2) in the abovementioned example. The obtained values are subsequently plugged into the expressions given in theorem 3.1 to find the steady-state values of the intermediates. That is, for matching rate constants between the core and an extension model, the solutions to the steady-state equations of the core model completely determine the solutions to the steady-state equations of the extension model.

The conservation laws, however, impose different constraints on the steady-state solutions for given conserved amounts. Specifically, by inserting (3.1) into (2.2), we obtain 3.3The steady states of the extension model solve (3.2) and (3.3), whereas they solve (3.2) and (2.3) in the core model. Equation (3.3) is nonlinear in the concentrations of the core species. Nonlinear terms in the conservation laws can cause the two models to have substantially different properties. This is reflected in the example in §4.

Importantly, if the system has no conservation law, then addition of intermediates cannot alter any property of the core model at steady state. This will be the case, for instance, when production and degradation of all core species in the model are explicitly modelled.

## 4. An example

The number of steady-state solutions for the core model and an extension model can differ substantially. For matching rate constants, the steady states of each system are found by intersecting the steady-state equations for the core species with the conservation laws of each of the systems. The number of points in this intersection might differ between extension models and the core model, depending on the form of the conservation laws.

We illustrate this using the two-site phosphorylation system in figure 1*a* and include dephosphorylation reactions
4.1In addition, we add the reactions
4.2The motivation for the addition is not biological but for illustrative reasons. It allows us to plot the steady-state equations in two dimensions. We will consider the positive steady states of the core model in figure 1*a* together with (4.1) and (4.2), and the extension model in figure 1*c* together with (4.1) and (4.2), and *y*_{1} = *S*_{0} + *E*, *y*_{2} = *S*_{1} + *E*, *y*_{3} = *S*_{2} + *E*. The added reactions are core reactions and do not involve intermediates. Because the substrate *S*_{2} is degraded (*S*_{2} *→* 0), the total amount of substrate is no longer conserved and there is only one conservation law, namely that for the kinase (compare (2.3)).

At steady state, the core and any extension model fulfil the relation
4.3for some constants *a*_{1}, *a*_{2} > 0 that depend on the rate constants of each model (see electronic supplementary material).

The relation is obtained from the steady-state equations of the core model alone and therefore must be fulfilled by all extension models for matching rate constants (theorem 3.2). One can show that the concentrations of [*S*_{1}] and [*S*_{2}] at steady state are uniquely determined by [*E*] and [*S*_{0}] (see electronic supplementary material).

For a given conserved amount for the kinase, the steady-state concentrations are determined by the common points of the graph of (4.3) and the curve for the conservation law. For the core model, this curve is , which is a vertical line in the ([*E*], [*S*_{0}])-plane. Because (4.3) is strictly decreasing in [*E*], it follows that there is a single steady state for any choice of (figure 2*a*).

Consider next the extension model corresponding to figure 1*c* (with the modifications introduced in (4.1) and (4.2)). For arbitrary fixed rate constants of the core model, we choose rate constants of the extension model that realize the rate constants of the core model. This can always be achieved for extension models with ‘dead-end’ complexes, like that of figure 1*c* (see electronic supplementary material). For , the concentration of the intermediate is [*Y*] = *a*_{3}[*E*][*S*_{0}] for some constant *a*_{3} > 0 that depends on the rate constants of the extension model. Consequently,
4.4If *a*_{3} = 0, then we obtain the core model. In the particular case, *a*_{3} varies independently of *a*_{1}, *a*_{2} and all values of *a*_{3} can be obtained when realizing the rate constants of the core model. Combining (4.3) and (4.4) yields a second-order polynomial in [*E*]:
4.5Hence, for fixed *a*_{1}, *a*_{2}, *a*_{3}, the polynomial can have zero, one or two positive solutions, depending on the value of *E*_{cons}^{C}. Figure 2 shows graphically the steady-state solutions for the core (figure 2*a*) and the extension (figure 2*b*,*c*) model as the intersection of the steady-state equation (4.3) and the curve for the conservation law for different values of and *a*_{3}. In figure 2, the curve for the steady-state equation (dashed) is given for *a*_{1} = 2 and *a*_{2} = 0.5 and is the same for the two models. For the core model, the conservation law curve is a vertical line (solid, figure 2*a*), which intersects the steady-state curve in precisely one point (figure 2*a*). For the extended model, the conservation law curve (solid, figure 2*b*,*c*) is the ratio in (4.4). Depending on the value of , the two curves intersect in zero, one or two points illustrating how the number of steady states vary with (figure 2*b*, with *a*_{3} = 2). The same conclusion is obtained by varying *a*_{3} while keeping fixed (figure 2*c*, with ).

In this particular case, we could find explicit expressions for the steady-state concentrations in terms of the conserved amounts and the rate constants. This is not always the case.

## 5. Number of steady states

In the example in §4, one can choose rate constants and conserved amounts such that the extension model does not have a positive steady state, even though the core model has a positive steady state for all choices of rate constants. However, it is easy to see that *a*_{3} can always be chosen so small that there is at least one positive solution for fixed *a*_{1}, *a*_{2} and . If *a*_{3} ≈ 0, then the contribution of [*E*][*S*_{0}] in (4.4) becomes insignificant, and the extension model is ‘similar’ to the core model. This is observed in figure 2*c* for small *a*_{3}, the curve for the conservation law is almost a vertical line.

Therefore, in the example, it is always possible to choose matching rate constants such that the number of steady states in the extension model is at least as big as the number of steady states in the core model, for corresponding conserved amounts. This observation holds generally. We now state the main result concerning the dynamical properties of extension models and the number of steady states.

### Theorem 5.1.

*If the core model has N non-degenerate*^{1} *positive steady states for some rate constants and conserved amounts, then any extension model that realizes the rate constants has at least N corresponding non-degenerate positive steady states for some rate constants and conserved amounts. Oppositely, if the extension model has at most one positive steady state for any rate constants and conserved amounts then the core model has at most one positive steady state for any matching rate constants and conserved amounts. The rate constants and conserved amounts can be chosen such that the correspondence preserves unstable steady states with at least one eigenvalue with non-zero real part and asymptotical stability for hyperbolic steady states*.

The proof essentially relies on the observation in the previous example that a certain parameter (*a*_{3} in the example) can be chosen so small that the extension model and the core model are almost identical at steady state. The relationship between a reaction network and a subnetwork has been studied previously, but in different contexts. For example in [23,24], where subnetworks are defined by (certain) subsets of reactions, or in [24], where subnetworks are defined by removing species from reactions. Characterizations similar to theorem 5.1 about the number of steady states hold in these situations.

In figure 2, the steady state in the extension model closest to the steady state in the core model (for the same conserved amount) inherits the stability properties of the steady state of the core model. In this case, it is asymptotically stable. However, we cannot conclude anything about the other steady state in the extension model from the core model alone.

## 6. Steady-state classes and canonical models

The observations made about the conservation laws and the steady-state equations (theorems 2.1, 3.1 and 3.2) suggest that it suffices to know what core complexes contribute to the conservation laws in order to compare the extension and core models at steady state. In figure 1*b*, the core complexes *S*_{0} + *E*, *S*_{1} + *E* contribute to the conservation laws for the kinase and substrate. Any other extension model, contributing the same core complexes to the conservation law, will result in equations for the steady states of the same form. Specifically, if two extension models contribute the same core complexes to the conservation laws and realize the same rate constants,^{2} then the two models are identical at steady state. In particular, we can apply theorem 5.1 to any of the two models.

Therefore, we can group extension models according to the core complexes that appear in the conservation laws. We say that two extension models belong to the same *steady-state class* if they share the same core complexes in the conservation laws. The complexes characterizing a steady-state class are called the *class complexes*. We can use theorem 3.1 to provide a graphical characterization of the classes: the core complexes that contain a species appearing in some conservation law are selected. If there exists a reaction from such a core complex to an intermediate, then the core complex is a class complex. The class of the core model is the class with no class complexes.

In figure 3, the graphical characterization is illustrated using the core model in figure 1*a*, written in simplified form. All species appear in some conservation law, and hence all core complexes can be class complexes. Consider the extension models in figure 1*b* and figure 1*c,d* with *y*_{1} = *S*_{0} + *E*, *y*_{2} = *S*_{1} + *E*, *y*_{3} = *S*_{2} + *E*. The extension model in figure 1*c* belongs to the steady-state class with class complex *y*_{1}, because there is only one path from a core complex to an intermediate: *y*_{1} → *Y*. Similarly, the extension models in figure 1*b,d* have class complexes *y*_{1}, *y*_{2}. We conclude that figure 1*b,d* are in the same class, whereas the models in figure 1*a,c* are in different classes and have different equations. In this case, figure 1*a* has always one steady state for any choice of conserved amounts and rate constants, whereas figure 1*b–d* can be multistationary (this is proved by direct computation of the steady states in the electronic supplementary material).

Because class complexes characterize the steady-state classes, there is a finite number of classes, at most 2^{K}, with *K* the number of core complexes (*K* = 3 in figure 3). The classes are naturally ordered by set inclusion: a class is smaller than another class if the latter contains the class complexes of the former. In particular, the steady-state class of the core model is smaller than any other class. Thus, the class of figure 1*a* is smaller than the classes of figure 1*b–d*, then classes of figure 1*b,d* are the same and the class of figure 1*c* is smaller than the class of figure 1*b*. The classes of the models in the first and the third box of figure 3 are not comparable as the first is {*y*_{1}} and the last is {*y*_{2}, *y*_{3}}.

All extension models in a steady-state class have common properties at steady state (subject to the requirement of realizability of rate constants). Thus, it is natural to select a representative for each class with a small number of intermediates and such that the behaviours of all models in the class are reflected in the behaviour of the representative. To each class, we construct a *canonical model* by adding a dead-end reaction, , for each class complex (see figure 3 for an example). Importantly, the steady-state equations for the canonical model are simpler than for any other extension model in the same class. It is shown in the electronic supplementary material that the parameter space of the canonical model is as large as possible. This leads to the following corollary to theorem 5.1.

### Corollary 6.1.

*If the canonical model of a steady-state class has a maximum of N steady states for any rate constants and conserved amounts, then all extension models in the class, or in any smaller class, have at most N steady states.*

In particular, if the largest canonical model (with a dead-end reaction added to all core complexes) is not multistationary, then no extension model, including the core model, can be multistationary. Likewise, if the smallest canonical model (the core model) is multistationary, then all extension models are multistationary. If there are no conservation laws, then there is only one steady-state class and any steady state in the core model corresponds precisely to a steady state in the extension model (assuming rate constants are realizable; theorems 3.1 and 3.2). Hence, either all extensions models (with realizable rate constants) and the core model are multistationary or none of them are. Further, if the core model cannot have multistationarity neither can an extension model, independently of the realizability of the rate constants.

Theorem 5.1 and corollary 6.1 provide assistance to the model builder. First of all the modeller can focus on the canonical models only. By screening the canonical models for the possibility of multistationarity, the modeller obtains a clear idea about the effects of intermediates. In figure 3, the steady-state class given by {*y*_{2}, *y*_{3}} does not have multiple steady states, hence the same holds for the classes {*y*_{2}}, {*y*_{3}} and the core model (theorem 5.1). Multistationarity in figure 3 (two first columns) is due to the nonlinearity introduced by [*y*_{1}] in the conservation laws, irrespective of the presence or absence of [*y*_{2}] and [*y*_{3}].

Our approach provides a simple graphical procedure to classify the extension models into a finite set of classes with common dynamical features, thereby elucidating the consequences of choosing a specific model. Figure 4 shows a decision diagram that guides the modeller through a number of possibilities. Each decision can be checked using various computational methods [25–28] or by manually solving the system (a task that simplifies due to the simple form of the canonical models).

## 7. Example: two-component systems

Table 1 shows a biological application of the decision diagram in figure 4. We consider three models of two-component systems of increasing complexity [29,30]. The basic mechanism consists of a sensor kinase that autophosphorylates SK⇋SK* (here asterisk indicates a phosphate group), the phosphate group is subsequently transferred to a response regulator RR and dephosphorylation of RR^{*} is catalysed by a phosphatase Ph. This model is considered in table 1 (model A). Models B and C in table 1 consist of the first model enriched with more mechanisms. In model B, SK has a bifunctional role and acts as a phosphatase, and likewise RR catalyses dephosphorylation of SK. Model C is an enrichment of model B with dephosphorylation of SK^{*} by a phosphatase T. Models B and C in table 1 are core models of the models considered in [29,30]. Models B and C are not extension models of model A, nor of each other. All models considered in table 1 have the total amount of kinase and the total amount of response regulator conserved.

We have applied the decision tree in figure 4 to each of the models. Models A and C are *robust* with respect to the choice of intermediates: model A cannot exhibit multistationarity for any choice of rate constants and model C exhibits multistationarity for some choice of rate constants, independently of how intermediates are included in the models. On the other hand, model B is sensitive to how intermediates are introduced. The core model is not multistationary, but inclusion of intermediates in some reaction paths introduces multistationarity. We conclude that modelling of this system needs to be performed carefully, as the qualitative conclusions that can be drawn from the model depend on the choice of intermediates.

Our analysis of the canonical models identifies the steady-state classes that can exhibit multistationarity and pinpoints the particular class complexes that introduce nonlinearity in the conservation laws. The analysis provides a simple overview of the effect of introducing intermediates in different reactions.

## 8. Discussion

Our work develops from the perspective of the model and clarifies the effects of intermediate species in biochemical modelling. Simplifications are always applied in model building but generally on a case to case basis, motivated by biological assumptions. One example is the QSSA, where equations of the form , together with some (but not all) conservation laws, are used to eliminate species [20,22]. This results in a *hybrid* model between our core and extension models. Our framework allows us to eliminate interme diate species generally and to compare core and extension models in a formal mathematical way. This comparison can be made independently of particular biological assumptions. An important insight is that model simplification and model choice must be pursued with great care as crucial dynamical properties might change radically by the inclusion of intermediates.

We remarked in the introduction that intermediates have been shown to affect steady-state properties of a system, such as the emergence of ultra-sensitivity [5,6]. It follows from our results that intermediates cannot change a model's properties at steady state if there are no conservation laws. In particular, if production and degradation of each species are explicitly modelled, then a model without intermediates is fully justified at steady state.

It has previously been noted that models that seem very similar can have different qualitative properties [31]. Our analysis is a step forward in quantifying the relationship between simple and complex models of the same system, and in using simple models to predict properties of complex systems. Our results can guide the modeller through the critical issue of choosing a model and in learning about model properties. As such, the results are useful for interpretation of experimental data and for designing synthetic systems. We envisage that our techniques can be extended to other models than those defined by intermediates and can provide further insights into the nature of biochemical and other types of modelling [17,32].

## Funding statement

This work was supported by the Lundbeck Foundation, the Leverhulme Foundation and the Danish Research Council. E.F. is supported by a postdoctoral grant ‘Beatriu de Pinós’ from the Generalitat de Catalunya and the project no. MTM2012-38122-C03-01 from the Spanish Ministerio de Economía y Competitividad.

## Acknowledgements

Part of this work was carried out while E.F. and C.W. visited Imperial College London in 2011. Neil Bristow is thanked for assistance. The anonymous reviewers are thanked for their constructive comments.

## Footnotes

- Received May 30, 2013.
- Accepted July 1, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.