Parameter inference to motivate asymptotic model reduction: an analysis of the gibberellin biosynthesis pathway

Developing eﬀective strategies to use models in conjunction with experimental data is essential to understand the dynamics of biological regulatory networks. In this study, we demonstrate how combining parameter estimation with asymptotic analysis can reveal the key features of a network and lead to simpliﬁed models that capture the observed network dynamics. Our approach involves ﬁtting the model to experimental data and using the Proﬁle Likelihood to identify small parameters and cases where model dynamics are insensitive to changing particular individual parameters. Such parameter diagnostics provide understanding of the dominant features of the model and motivate asymptotic model reductions to derive simpler models in terms of identiﬁable parameter groupings. We focus on the particular example of biosynthesis of the plant hormone gibberellin (GA), which controls plant growth and has been mutated in many current crop varieties. This pathway comprises two parallel series of enzyme-substrate reactions, which have previously been modelled using the law of mass action [23]. Considering the GA20ox-mediated steps, we analyse the identiﬁability of the model parameters using published experimental data; the analysis reveals the ratio between enzyme and GA levels to be small and motivates us to perform a quasi-steady state analysis to derive a reduced model. Fitting the parameters in the reduced model reveals additional features of the pathway and motivates further asymptotic analysis which produces a hierarchy of reduced models. Calculating the Akaike information criterion and parameter conﬁdence intervals enables us to select a parsimonious model with identiﬁable parameters. As well as demonstrating the beneﬁts of combining parameter estimation and asymptotic analysis, the analysis shows how GA biosynthesis is limited by the ﬁnal GA20ox-mediated steps in the pathway and generates a simple mathematical description of this part of the GA biosynthesis pathway.


Introduction
Ordinary differential equation (ODE) models of biological regulatory networks usually contain many unknown parameters. Such parameters are typically very challenging to measure experimentally and are often estimated indirectly by calibrating the model to experimental data such as a time course of measurements under an experimental perturbation. Such model calibration deepens understanding of the network dynamics and enables the model to make more accurate predictions (predicting, for instance, the dynamics of unseen components or the behaviour of the network under alternative perturbations). Furthermore, the calibration enables inference to be made about the parameter values themselves, which are often of scientific interest.
The question of whether or not parameters can be estimated well, or at all, from experimental data is called "identifiability". There are many good reviews of approaches for determining identifiability of ODE models (see e.g. [22]) but relatively much less on how models with non-identifiable parameters can be reduced to obtain iden-tifiable models, which is the focus of this paper.
A characteristic common to many ODE models is that model dynamics and data agree for parameter values that range over many orders of magnitude, and hence the parameters are not precisely identifiable from data. This phenomena is often referred to as "sloppiness" [4,14] (which we will define more precisely below). In some circumstances, sloppiness can be mitigated with careful experimental design, i.e., the choice over what experiments to perform, including what to measure and when [30]. There are also circumstances in which sloppiness causes no issues; for example, even though parameter estimates may have very high variance, predictions for particular variables can still have low variance [4,14].
In contrast to these approaches, here we consider how to reduce a sloppy model to a simpler one with identifiable parameters. We show how analysing parameter sensitivity provides understanding of the model dynamics that enables us to identify suitable asymptotic model reductions. Although the full model encapsulates the key mechanisms thought to generate the data, the reduced model is an approximation that describes the experimental data nearly as well. The reduced model reveals the level of complexity required to capture the dominant dynamics on a particular time scale, and suggests how these dynamics depend on groupings of the original model parameters. Deriving the reduced model retains mechanistic interpretation and physical meaning of the parameter groupings, and fitting the reduced model reveals how well these parameter groupings can be identified from the data. Deriving reduced models is helpful not only for parameter inference but also for computational reasons: complex multiscale models are often created by coupling submodels for different components of a large system (e.g. [7,3]) so having smaller model components helps to reduce the overall complexity and computational cost.
The particular application we consider is the biosynthesis of the plant hormone gibberellin (GA), which regulates organ growth, flower development and seed germination, and plays a key role in adaptation to environmental stress [8,11]. GA biosynthesis has been shown to be predominantly regulated at the later steps in the pathway whereby GA 12 is converted to active forms GA 4 and GA 1 [12]. Here, the GA biosynthesis pathway comprises two parallel series of enzyme-substrate reactions [15,17,33] (shown in Fig. 1): GA 12 either undergoes a series of reactions catalysed by members of the GA20ox family, to produce GA 9 which is then converted to the active GA 4 via the GA3ox enzymes, or GA 12 is converted to GA 53 (via GA13ox) and then undergoes a series of GA20ox-mediated steps to produce GA 20 which is converted to the active GA 1 by GA3ox. The active GAs and their precursors can also be degraded via the GA2ox family of enzymes. Thus, the level of bioactive GA depends on the balance between synthesis and degradation, which are controlled via the levels of GA20ox, GA3ox and GA2ox enzymes; these enzymes exhibit distinct spatial and temporal patterns and are regulated by environmental signals and other plant hormones [15,17,18,25,33].
The GA pathway has played a major role in improving crop productivity: during the Green Revolution in the 1960s, high-yielding semi-dwarf cereal crop varieties were developed, many of which were later found to have genetic mutations in the GA pathway [26]. Understanding the mechanisms that control GA levels is a prime target for further increasing yields [8,16].
In this study, we focus on the GA20ox-mediated steps of the GA biosynthesis pathway and describe how parameter estimation, parameter sensitivity analysis and asymptotic analysis can generate a simpler description of this biological network, providing understanding of the underlying dynamics and identifying the parameter groupings that govern the network dynamics. We first describe our methods for a general case ( §2), describing fitting criteria, derivation of confidence intervals and how the profile likelihood can be used to motivate asymptotic model reduction. We then apply our methods to analyse the dynamics of the GA20ox-mediated steps in the GA biosynthesis pathway ( §3). Beginning with the model suggested by Middleton et al. [23], based on the law of mass action, we develop a series of reduced models, assessing each model by comparing simulation results with the experimental data from Appleford et al. [2]. We demonstrate how our methods motivate the derivation of a simple description of the observed dynamics and reveal the key network interactions. Our conclusions are discussed in §4.

Models and fitting criteria
The dynamics of a biological network are often modelled by a system of ODEs of the form where y = y(t; θ) = (y 1 (t; θ), . . . , y s (t; θ)) is a vector containing variables representing concentrations of the various interacting species included in the model; θ = (θ 1 , . . . , θ p ) is a vector of models parameters; and y 0 is the initial value of y at time t = 0. We suppose that only a subset of the y 1 , . . . , y s -here taken without loss of generality to be the first m ≤ s -are measurable experimentally, and that these are measured only at a finite number of time points and subject to measurement error. Supposing that variable y j , say, is measured at the n j time points t j1 , . . . , t jnj (not necessarily the same for different j) and that y Data ji denotes the experimental measurement of variable y j at time t ji , then the full set of experimental data, comprising N = n 1 + · · · + n m data points, is y Data ji ; i = 1, . . . , n j , j = 1, . . . , m .
Fitting a model (1) to experimental data (2) involves determining θ (and possible some elements of y 0 if they themselves are unknown) so that the model solution is close, by some criterion, to the data. One approach is to minimise the sum-of-squares criterion, to findθ = argmin S(θ), which is the "least-squares estimate" of θ. Another, more statistically minded, approach is to define an error model relating the deterministic model (1) to the data (2), construct the "likelihood function", L(θ) (see e.g. [9]) then estimate θ by maximising L(θ) (or equivalently the loglikelihood function, (θ) = log L(θ)). For example, under a Gaussian error model where the errors ε ij ∼ N (0, σ 2 ) are independent Gaussian random variables (i.e. Normally distributed) with mean 0 and variance σ 2 , the log-likelihood function is where N is the number of data points and S(θ) is given by (3). For any σ 2 , (θ) is maximised with respect to θ by minimising S(θ). In other words, under a Gaussian error model, (5), which we adopt throughout this paper, the "maximum-likelihood estimate",θ = argmax (θ) is identical to the least-squares estimate of θ (4). The advantage of the likelihood framework, however, is that it provides distributional results for making statistical inference about θ, as described below.

Confidence intervals for parameter values
The asymptotic distribution (i.e., as N → ∞) ofθ, under mild conditions that commonly hold, iŝ where N p (·, ·) denotes the multivariate Gaussian distribution, and I(θ) = −EH is the so-called Fisher information matrix, where E denotes expectation and is the Hessian matrix of with respect to θ [29]. Result (7) can be used to construct individual or joint confidence regions for the components of θ. For example, from (7), by basic properties of the multivariate Gaussian distribution [21],θ j is distributed aŝ then, approximating I(θ) with −H(θ), we can derive the bounds of a 100(1 − α)% confidence interval aŝ where χ 2 1 (α) is the α upper quantile of the χ 2 1 distibution. The term "sloppiness", introduced above, is now sometimes used very loosely to mean "difficulty in estimating parameters", though the term was originally introduced with the more particular meaning that the spectrum of eigenvalues of H span many orders of magnitude, and include some very small values. Small eigenvalues make H ill-conditioned, inflating elements of H −1 and hence inflating the confidence intervals for the individual parameters. We show later in Fig. 7 the eigenvalue spectra of models of different complexity, and in Tables 2 and 3 examples of confidence intervals computed using (9). For details of the numerical calculations of H in this paper, see Appendix A.
A different way to investigate the information about particular parameters contained within the data is by the "profile likelihood" [28,32]. This approach explores parameter space in the neighbourhood ofθ by varying a single parameter and fitting the model with respect to the remaining parameters. As described below, this provides insight into the shape of the likelihood in the neighbourhood of the maximum. To calculate the profile likelihood, consider a particular parameter of interest, θ j , denote the other parameters, θ −j , and write the log-likelihood as (θ j , θ −j ). The log-profile-likelihood is defined as whereθ −j = argmax θ−j (θ j , θ −j ), i.e., (10) is the log likelihood maximised over all the parameters except the particular one of interest. By Wilks' theorem [31], under a null hypothesis that θ j takes a particular value, θ 0 j , say, then asymptotically which can be rearranged to give the confidence interval The confidence interval can be visualised by marking the χ 2 1 (α) threshold on a plot of −2 P (θ j ) versus θ j ; see for example, Fig. 10.
An intuitive view of the Gaussian-and profile-likelihoodbased confidence intervals, (9,12), is that they both involve investigating the curvature of (θ) in the vicinity of θ =θ: flatness indicates that the data are not informative about θ, whereas high curvature indicates that the data are very informative about θ. The Gaussian-based confidence intervals amount to approximating (θ) by a quadratic with its curvature matched at the maximum,θ. This entails an assumption of symmetry that may not be warranted for the model and data (see, e.g., Fig 11). The profile-likelihood confidence intervals involve no such assumption and are often more accurate than the Gaussianbased confidence intervals for nonlinear models and data sets of typical size [28,32]. Another advantage of profilelikelihood confidence intervals is that they are invariant to the choice of parameterisation, e.g. whether or not the parameters are defined on a log scale (as in this paper), whereas Gaussian-based ones are not.
A point sometimes not highlighted is that the asymptotic results (7,11) and hence confidence intervals derived from them, (9,12), require that θ is a point in the interior of parameter space rather than a point on the boundary. This is not always the case: for example, some of the parameters that are non-negative on physical grounds may be exactly zero. In Model II later (θ) is maximised when one or more of the parameters approaches zero, in which case we do not compute confidence intervals. Another point we note is that in the derivations of the confidence intervals above we have assumed that σ 2 is known, whereas here we substitute σ 2 with its estimateσ 2 = N −1 S(θ), where S(θ) is computed from the fitted model. This has little effect on confidence intervals and saves the expressions for the confidence intervals from becoming more complicated. Details of how to compute the confidence intervals without making this substitution are given in [29, p196].

Profile likelihood motivating model reduction
As we demonstrate in our case study below, the presence of non-identifiable parameters suggests that the observed dynamics can be captured by a simpler model. To derive the simpler model, the profile likelihood can be used to suggest an appropriate asymptotic model reduction. For example, if the confidence interval is unbounded in one direction, we take the corresponding limit in which θ i → 0 or θ i → ∞. In some cases, the model dynamics may depend on θ j only in a ratio with other parameters, in which case plotting how the other fitted parameters,θ −j , vary as θ j is varied suggests an appropriate model reduction. Taking an asymptotic limit may also be appropriate for an identifiable parameter if the estimated value and profile likelihood show this parameter to be large or small. Having derived a reduced model, we compute an information criterion (AICc, see Appendix B) to select amongst the candidate models. A schematic of this model reduction approach is shown in Fig. 2.

Modelling the GA20ox-mediated steps in the GA biosynthesis pathway
The GA biosynthesis pathway features two parallel series of oxidation steps mediated by the GA20ox enzymes ( Fig. 1). Each of these steps involves the GA species first binding to GA20ox with a reversible reaction, and the resulting complex then dissociating into the next GA species in the pathway and GA20ox [23] (Fig. 3a). Thus, the con-version of GA 53 to GA 20 (Fig. 3b) involves the reactions where k dj , k aj and k mj denote rate constants for j = 53, 44, 19, square brackets denote concentrations and asterisks denote dimensional quantities. Equivalent reactions govern the conversion of GA 12 to GA 9 (Fig. 3c); however, to maintain conciseness, we give reactions and equations for only the GA 53 pathway throughout.

Middleton et al's model: Model I
To describe the reactions (13), we initially adopt the model by Middleton et al. [23] who modelled the conversion of GA 12 to GA 9 using the law of mass action. Thus, reactions (13) can be described by a system of coupled  ordinary differential equations (ODEs): where t * denotes time.
To estimate the parameters in this model, we follow [23] in comparing the model dynamics to data from Appleford et al. [2] who performed in vitro studies to investigate the dynamics of the conversion of GA 53 to GA 9 and GA 12 to GA 20 . Considering the GA 53 conversion, for ex-ample, their solution initially contained fixed concentrations of GA 53 and GA20ox1, which we denote by s 0 and e 0 respectively, and the concentrations of GA 53 , GA 44 , GA 19 and GA 20 were measured at subsequent time points (see data points in Fig. 4). Thus, the data corresponds to the initial conditions [ 19 .GA20ox] * = 0. In the experiments [2], the initial concentration s 0 is given, which following [23] we take to be known initial conditions. The initial GA20ox concentration e 0 is an unknown and is a parameter to be estimated.
Fitting the model parameters to the data (as described in §2.1 and Appendix A), we see that equations (14) (hereafter referred to as Model I), can represent the observed dynamics for the GA 53 and GA 12 pathways; see blue lines in Fig. 4 and 5 respectively (consistent with the results in [23] for the GA 12 pathway). The profile likelihood for each parameter (Fig. 6) shows that the fitted model remains in agreement with the data for wide ranges of parameter values. This conclusion is consistent with the eigenvalue spectrum of the Hessian varying over several orders of magnitude (Fig. 7), and hence the model being sloppy. The model's sloppiness motivates us to seek a simpler model for the observed dynamics. The profile likelihood suggests that to fit the data the estimated initial GA20ox concentration, e 0 , needs to be small compared to the prescribed initial GA 53 concentration, s 0 . Thus, motivated by the profile likelihoods and previous studies of enzymesubstrate reactions, we will consider the limit in which e 0 s 0 .

Model reduction, Model II
Applying the law of mass action is a standard way of physically describing enzyme-substrate reactions; however, the resulting Model I is complex with the model dynamics depending on 10 parameter values. Given the profile likelihood revealed that Model I only fits the data if e 0 /s 0 is small, we now consider a model reduction by taking the limit e 0 /s 0 → 0. This quasi-steady-state assumption is the basis of the Briggs-Haldine derivation of the Michaelis-Menten equation [24], and, for a single substrate-enzyme reaction, has been shown to generate a reduced model in which parameters are identifiable using standard time course data [6]. In the limit as e 0 /s 0 → 0, the complexes form rapidly and the complex concentrations are in equilibrium after an initial transient time scale.

Nondimensionalisation
To derive a reduced model, we first non-dimensionalise the governing equations. We consider the GA concentrations relative to the initial GA 53 concentration, s 0 , the complex and enzyme concentrations relative to the initial enzyme concentration, e 0 and time relative to the rate at For these data, the known initial concentration of GA 53 is s 0 = 5.56 µM. We note that there is no noticeable difference between the predicted Model II and Model III dynamics and that in this case, Model IV does not capture the dynamics well (as is reflected in the S(θ) and AICc values in Table 1). In the upper panels, the predicted GA 53 , GA 44 , GA 19 and GA 20 for Models II-IV are redimensionalised via (15) to show comparison with the experimental data.
We introduce the following dimensionless parameter groupings: On nondimensionalising in this fashion, the governing equations, (14), become

Model reduction
Having non-dimensionalised, we reduce the model by taking the limit e 0 /s 0 → 0. Considering t * = O(1), equations (17e-g) give expressions for the complex concentrations in the limit as → 0 where (i.e. κ j = k aj s 0 /(k dj + k mj ) the ratio between the typical complex association rate and the complex dissociation rate). Summing (17e-h) shows that the total enzyme concentration is conserved; thus the initial conditions, (18b-c), require that the enzyme and complex concentrations sum to one for all times, and from (19), we obtain .
where parameter groupings These equations and expressions depend on 6 parameter groupings (20,23). The reduced model cannot satisfy initial conditions (18b-c). As in the Michaelis-Menten approximation of a single enzyme-substrate reaction, the complex concentrations evolve to their quasisteady state on an initial time scale 0 ≤ t * 1 [24]; however, to remain focussed on data-fitting, we consider the derivation of this inner solution to be beyond the scope of this study.
Fitting Model II to the experimental data [2] (as described in Appendix A), we see that the reduced model captures well the observed dynamics ( Fig. 4 and 5), suggesting that the complex concentrations are indeed quasisteady on the time scale of the experiment. Since Model II is derived as an asymptotic reduction of Model I, it fits the data slightly less well than Model I for both pathways in the sense of having a larger residual sum-of-squares, S(θ), but it involves four fewer parameters and by AICc Model II is favoured over Model I (see Table 1 and Appendix B).
The profile likelihoods show different behaviour for the two pathways (Figs. 8 and 9). For the GA 53 pathway, the profile likelihood indicates that the model dynamics fit the data well when κ 53 and κ 44 are small; however, the confidence intervals of the remaining four parameters are finite (and do not include zero) (Fig. 8). For the GA 12 pathway, we similarly see that κ 12 and κ 15 are small, but in this case, the other parameters have unbounded confidence intervals and hence any of these may be very large (Fig. 9). For both pathways, the eigenvalue spectra of the Hessian is sloppy (Fig. 7). Given the similar features in the profile likelihoods of the two pathways, we consider a model reduction by taking κ 53 and κ 44 (or equivalently κ 12 and κ 15 ) to be small. We note that the profile likelihoods (Figs. 8 and 9) also show that the fitted values of the other parameters do not change as κ 53 , κ 44 , κ 12 and κ 15 varied from small values, and thus, it appears that these four parameters have a negligible influence on the model dynamics on this time scale.

A further model reduction: Model III
The profile likelihoods ( Fig. 8 and 9) show that the model fits the data equally well for any small values of the parameters κ 53 and κ 44 for the GA 53 pathway and κ 12 which are coupled to an expression for the enzyme concentration and expressions (19) for the complex concentrations, and subject to initial conditions (24). In this limit, the dissociation rates of complexes GA 53 -GA20ox and GA 44 -GA20ox are much larger than their association rates so that their concentrations are asymptotically small. Fitting the parameters for Model III, we find no noticeable difference between the predicted Model III and Model II dynamics (see Fig. 4 and 5); the fitted parameters and residual sum-of-squares, S(θ), are similar to those of Model II; however, since Model III has fewer parameters than Model II, Model III has a lower AICc (Table 1). For both pathways, the range of eigenvalues of the Hessian is much smaller for Model III than Model II, showing that Model III is less sloppy that Model II (Fig. 7). The profile likelihoods show that the parameters are identifiable for the GA 53 pathway and we find reasonable agreement between the Gaussian-and profile-likelihood-based confidence intervals in this case (Fig. 10, Table 2). In contrast, the profile-likelihood based confidence intervals are unbounded for the GA 12 pathway (Fig. 11). Thus, the profile likelihoods suggest that dynamics of the GA 53 pathway appear to be parsimonously described by Model III, whereas those of the GA 12 pathway may be captured by a simpler model. The profile likelihoods for the GA 12 pathway show that any of the four parameters may be large, and there appears to be a linear relationships between the parameters along each profile likelihood (Fig  11), which suggests that the model dynamics depend on ratios between certain parameter values.

A final model reduction: Model IV
Motivated by the profile likelihoods of the GA 12 pathway, we consider the limit in which κ 19 → ∞, and introduceλ j = λ j /κ 19 for j = 53, 44, 19 which we take to be and We now consider the limit as δ → 0. In this case, the leading-order equations depend on the magnitude of [GA 19 ] compared to δ, therefore to derive the governing equations, we consider a series of distinct timescales in which different terms are present at leading order (see [10,19] for examples where similar approaches have been used, and [20] for a textbook description). As detailed below, five timescales emerge from the asymptotic analysis; to validate this analysis, Figs. 12 and 13 compare the asymptotic solutions to numerical simulations of Model III, (27), with δ = 0.01.

3.4.1.
Matched asymptotic analysis to evaluate the leading-order dynamics in the limit δ → 0 Timescale 1, t * = O(δ 2 ): At the beginning of the experiment, only GA 53 and GA20ox are present. Therefore at early times [GA 19 ] δ, the amount of enzyme bound with [GA 19 ] is asymptotically small and the concentration of free enzyme [GA20ox] = 1 (see (28)). We rescale (using daggers to denote the rescaled variables) and find that the concentrations evolve on a fast time scale, T , such that t * = δ 2 T . At leading order, equations (27) (satisfying initial conditions (24)). As shown in Fig. 12a, these analytical expressions agree with the dynamics of Model III, (27), at early times.
Timescale 2, t * = O(δ 3/2 ): As the GA 19 concentration increases, we move to a different asymptotic regime in which [GA 19 ] = O(δ) and significant amounts of both free and GA 19 -bound enzyme are present. To ensure matching between consecutive timescales, we require the long-time behaviour of each variable to match to the short-term dynamics on the subsequent timescale. Thus, considering the long-time dynamics of Timescale 1, if T is scaled with δ α as T → ∞ (for unknown α) then from (31), we find that with δ 3α . We find a suitable balance setting α = −1/2. Thus, the rescaled time,t, is t * = δ 3/2t and the relevant scalings are The resulting rescaled problem reads  20 ] ‡ = 0, numerical solutions reveal that these equations, (33), faithfully reproduce the dynamics of Model III at early times (Fig. 12b).
Given we do not have an analytical solution on this time scale, we elucidate the long-time behaviour of the dynamics by supposing each variable takes the form αt η ast → ∞ (where α and η are unknown constants for each variable); substituting into the governing equations (33), we find that These long-time dynamics are in good agreement with numerical solutions of the governing equations on this time scale, (33) (Fig. 12b).   (Fig. 12c).
Matching to the subsequent time scale requires the dynamics in the limit t * → 1/λ 19 . Supposing each variable takes the form α(1/λ 19 − t * ) η , we find suitable balances in the governing equations, (35), with These analytical solutions accurately represent the Timescale 3 dynamics in the limit as t * → 1/λ 19 (Fig. 13).        (Fig. 13), the leading-order solution for Timescale 5 is simply given by

Fitting Model IV
The dominant (t * = O(1)) dynamics are given by the asymptotic solutions in Timescales 3 and 5, which can be solved subject to prescribing an initial non-zero concentration of GA 19 . For a given parameter set, we simulate the Timescale 3 equations for t * ∈ [0, 1/λ 19 ] (35), prescribing initial conditions using the long time dynamics of Timescale 2 (34), applied at t * = 10 −4 , and then prescribe the Timescale 5 equations for t > 1/λ 19 (40). Simulating Model IV in this way, we fit the three model parameters to the experimental data (Figs. 4 and 5). The profile likelihoods (Fig. 14) show that for Model IV each parameter is identifiable using the data, with narrow Gaussian and profile-likelihood-based confidence intervals that are in good agreement ( Table 3).
As expected given the parameter estimates for Model III, the value of S(θ) and the AICc (Table 1) show that the Model IV approximation (κ 19 → ∞) captures the observed dynamics of the GA 12 pathway, but is not appropriate for the GA 53 pathway. Thus, our study recommends that Model IV is the appropriate choice for the GA 12 pathway, whereas one would choose Model III for the GA 53 pathway.

Discussion
We have described how parameter estimation can motivate asymptotic reduction of a network model and lead to a simpler model that represents well the observed dynamics. We demonstrated the use of these methods to analyse a recently published model of GA biosynthesis [23]. When fitting sloppy systems biology models to experimental data, calculating the profile likelihood provides helpful insight into the nature of parameter space in the neighbourhood of the best estimated parameter set. The profile likelihood reveals when taking a parameter to be increasingly small or large is consistent with the model fitting well the observed dynamics, or where the model dynamics are insensitive to particular parameter values. This, together with plots of how the other fitted parameters change as the prescribed one is varied, helps identify an appropriate asymptotic model reduction. The simpler model can then be fitted to the data to assess whether it is an appropriate approximation. As asymptotically reduced models depend on groupings of the original model parameters, the reduced model has fewer parameters and is a more parsimonious representation of the observed dynamics. Asymptotic model reduction often involves careful consideration of distinct time scales (for instance, as in the reduction of Model III to Model IV), and therefore automating the model-reduction step is unlikely to be possible in many cases.
The profile likelihood also enables calculation of confidence intervals for each parameter. These tend to be preferable in practise compared with confidence intervals based on the approximate Gaussian distribution ofθ, especially so for strongly non-linear models where the Gaussian approximation is poor (see e.g. Fig. 11). However, we have found that in the final reduced models, where parameters are identifiable, the two types of confidence interval do agree quite closely (e.g. Fig. 10).
We have focussed here on using the profile likelihood for parameter inference and model reduction. Other approaches for considering dependence between model pa-rameters with a view to model reduction include investigating the eigenvectors of H that correspond to the smallest eigenvalues [22], or adopting a Bayesian paradigm [13] and investigating structure amongst posterior parameter samples. In our view, though, the profile likelihood offers the most direct way to understand the effect of varying the value of individual parameters, and hence identifying good candidate model reductions.
Our study focussed on analysing a model of the GA biosynthesis pathway. Given the importance of GA for plant development, a GA pathway model is likely to feature in future plant models and the existing GA pathway model [23] has already featured in a multiscale model of plant growth regulation [3] and a study on hormone crosstalk [1]. In assessing the required level of complexity for capturing the dynamics of GA biosynthesis, our analysis will simplify future plant models. We recommend that Model III with the parameter estimates given in Table 2 could be used to simulate the conversion of GA 53 to GA 20 whereas Model IV with the parameter estimates given in Table 3 captures the conversion of GA 12 to GA 9 .
For both pathways, the analysis revealed the typical enzyme concentrations to be significantly smaller that the typical GA concentrations; as a result, the dynamics can be captured by Model II in which complex and enzyme concentrations equilibriate rapidly and are at a quasisteady state after an initial transient time scale. The concentrations of complexes or free enzyme can be approximated as simply a function of the GA concentrations (19,21). Fitting the parameter groupings under this Model II approximation further revealed that for GA 53 and GA 44 (and equivalently GA 12 and GA 15 ) the ratios between the typical complex association and dissociation rates (κ 53 and κ 44 respectively) are small, resulting in small GA 53 -GA20ox and GA 44 -GA20ox complex concentrations. We note that these parameter groupings, κ 53 and κ 44 , arose through the first asymptotic model reduction to Model II and therefore the second asymptotic limit is only revealed once Model II is fitted to the data. We find a further approximation is applicable to the GA 12 pathway; here, fitting Model III to the data showed that binding of GA 24 to GA20ox to form GA 24 -GA20ox occurs much more rapidly than GA 24 -GA20ox dissociation, resulting in the majority of the enzyme being bound in a complex with GA 24 . Thus, the modelling suggests that dissociation of GA 24 -GA20ox is a major limiting step in the production of bioactive GA 4 .