EGFR Signaling

We begin with an illustrative example; consider a model of EGFR signaling due to Brown et al. [6]. The model contains 48 unknown parameters (reaction rates and Michaelis constants) that were originally fit to 63 data points [6]. These data gave limited measurements of activity levels a few proteins after EGF and NGF stimulation in conjunction with a few knockout perturbations. Notably, this model gave a reasonable fit to experimental (not simulated) data. Furthermore, the model also made accurate, falsifiable predictions for the system behavior under novel perturbation experiments that were subsequently verified. From this, we conclude that the model likely contained all of relevant biology and chemistry necessary for giving a mechanistic explanation of these observations. However, the parameters were largely unconstrained when fit to this data, with relative uncertainties ranging from a factor of 50 up to a factor of a million.

In order to constrain the parameter estimates, Apgar et al. [13] proposed a set of five experimental conditions specifically for the model of Brown et al. These experiments were selected from a candidate pool of more than 160,000 possible experiments that included EGF and NGF stimulations at various levels and in various combinations, as well as potential knockout and over-expression perturbations. The five optimal experiments were chosen according to a “greedy” algorithm to maximize the smallest eigenvalue of the FIM. Under these experimental conditions, the authors expected to estimate all parameters to better than 10% accuracy. The experiments were simulated and not physically performed.

Rather than actually perform the five experiments proposed by Apgar et al., we create a second model to act as a surrogate for reality. The EGFR system is well studied and a model that reflects our current understanding of the system [42] would contain far more than the 48 parameters of Brown et al. As a reasonable next step along the ladder of realism, we replace the Michaelis-Menten approximations of the Brown et al. model by the mechanistic model from which this approximation was derived. The mechanistic model of (Henri-)Michaelis-Menten is two step enzyme-catalyzed reaction obeying mass-action kinetics: E + S ⇌ ES → E + P. This change introduces several new parameters (from 48 to 70) as well as several new chemical species corresponding to the intermediate enzyme-substrate complexes.

Note that there are a large class of mechanisms that could result in the same approximate Michaelis-Menten equation, so implementing the mechanistic model described above also makes several other simplifying assumptions. For example, a mechanistic model could also be written as E + S ⇌ ES ⇌ EP → E + P, i.e., with an isomerisation step for the enzyme-substrate complex. Indeed, there is a hierarchy of refining approximations one could make to this model. Our choice represents what is likely the simplest next step in refining the mechanistic description of Brown et al. For brevity, we refer to the original model of Brown et al. that implements the Michaelis-Menten approximation as the approximate model. We refer to the model of the Michaelis-Menten mechanism with mass-action kinetics as the mechanistic model.

Using the model and parameter values of Brown et al., we simulate the experimental conditions from reference [6] and add random noise to generate an initial data set characteristic of that in reference [7]. We then fit the 70 parameter, mechanistic EGFR model to this initial data. The resulting fit is both sloppy and unidentifiable as can be seen in Fig 4 (second column). Notably, the 70 parameter model has 22 more eigenvalues than the 48 parameter model, but there is not a clear separation between the largest 48 and the smallest 22 eigenvalues. Furthermore, it is not possible to equate the largest 48 eigenvalues with the 48 parameters of the approximate model and the 22 smallest eigenvalues with the new parameter combinations introduced by the more detailed kinetics.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. FIM for the four EGFR models.

Both the approximate Michaelis-Menten kinetics and mechanistic mass-action kinetics are unidentifiable when fit to the data in reference [7]. Although the optimal experiments in reference [13] lead to an identifiable (but still sloppy) model for the approximate Michaelis-Menten kinetics, the mechanistic mass-action kinetics remain unidentifiable. Furthermore, the FIM of the mass-action model suggests that a minimal model should include at least 60 parameters to explain the expanded observations, i.e., the manifold has approximately 60 widths larger than the experimental noise. The approximate Michaelis-Menten kinetics do not contain all of the relevant physics. The red dashed line corresonds to a relative standard error of 1/e in the inferred parameters.

https://doi.org/10.1371/journal.pcbi.1005227.g004

Next, we simulate artificial data for each of the experiments suggested by Apgar et al. using the 70 parameter mechanistic model and parameter values estimated from the first fit. We add random noise to these simulations to create a second data set. This second data set, having come from the more complicated model, acts as a surrogate for real experimental data. We then fit the 48 parameter, approximate model to the second artificial data.

The use of Michaelis-Menten kinetics in models of protein networks is somewhat controversial. In many practical cases, such as the original experiments of Brown et al., the model seems to work. However, the actual condition that must be satified is shown rigorously in reference [43] to be (3) where [E] and [S] are the enzyme and substrate concentrations respecively and K_M = (k_r + k_cat)/k_f is the Michaelis constant. From this condition we see that it is not strictly necessary for [E] ≪ [S] as is often (incorrectly) asserted. However, this condition is derived for a single enzyme-substrate reaction in isolation and, as we will see below, is generally not a sufficient condition for the Michaelis-Menten approximation to hold in a network context. Nevertheless, Eq (3) suggests that the Michaelis-Menten model could be valid provided K_M is very large, even if E ∼ S.

In choosing parameters for the mechanistic model, we therefore choose parameters such that the combination K_M = (k_r + k_cat)/k_f is very large. This is always possible because the model is unidentifiable when fit to the initial data. We were able to choose parameter values for all reaction rates such that the Michaelis-Menten approximation would not give errors larger than 10% for any of the reactions given the conditions of the network. This means that there is no a-priori reason to think that the approximate model of Brown et al. would be a poor approximation to the mechanistic model for our parameter values.

Even when enforcing the constraint that K_M is large, we find that the approximate model cannot give a reasonable fit to the data generated by the mechanistic model for the Apgar experiments. This is illustrated in Fig 5 where the systematic errors in the fit are clearly much larger than random noise.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Fit of approximate Michaelis-Menten kinetics to mechanistic mass-action data.

Because the approximate model does not contain all of the relevant mechanisms for the expanded observations, it cannot give a reasonable fit (i.e., within the expected variance of the experimental noise) to all of the experiments simultaneously. We see here that several time series are fit quite badly which could guide a modeler in identifying the missing relevant mechanisms.

https://doi.org/10.1371/journal.pcbi.1005227.g005

Because the error in the fit is large, we conclude that approximate Michaelis-Menten kinetics do not contain the relevant mechanisms to explain the observations under the expanded experimental conditions. This can be seen from the FIM eigenvalues in Fig 4. In particular, observe that the FIM for the mechanistic mass-action model under the expanded experimental condition contains approximately 60 eigenvalues larger than 1, so that there are about 60 directions on the model manifold with widths larger than the experimental noise. This indicates that a minimal model would require about 60 parameters to fit this data; the 48 parameter model of Brown et al. is clearly insufficient. Furthermore, the parameters of the mechanistic mass-action model are unidentifiable for the experimental conditions of Apgar et al. If the approximate model were replaced by the mechanistic model, it would require another round of experimental design in order to find accurate parameter estimates.

Because the fit of the mechanistic mass action kinetics to the original experiment is unidentifiable, there are many possible parameter values that would give equivalent fits. We generated an ensemble of parameter values for the mechanistic model consistent with the original experiments. When these ensembles generate data for the optimal experiments, we consistently find that the best fits for the approximate model have large errors similar to that in Fig 5.

The discrepancy between the data and the fit in Fig 5 is reminiscent of over-fitting. This is not the case however. Over-fitting occurs when the fit to a training data set is very good (i.e., too good) so that the predictions on a test set suffer as a result. In this language, the fit in Fig 5 is the training set (not the test set), so this is not an instance of over-fitting. Rather it is a demonstrationg that the model cannot fit the data.

Quantifying Model Error in Sloppy Systems

Motivated by the example of optimal experimental design in EGFR signaling described above, we now propose a simple method to quantify model descrepancy when fitting data to approximate models.

We assume that parameters are estimated by least squares regression, although many of our methods and results generalize. Least squares regression is justified by the assumption that the data are generated from a model with additive Gaussian noise: (4) where d_i is the i-th data point, y_i is the model prediction for the i-th data point, θ is a vector of parameters, ξ_i is a random variable with zero mean and unit variance, and σ_i is the scale of the noise. The random variable may account for any intrinsic stochasticity in the system, e.g., thermal fluctuations in particle number or inconsistencies in experimental measurements. It may also accommodate systematic errors such as model approximations, i.e., mechanisms that have been left out of the model. The field of uncertainty quantification has begun to explore methods to account for the inadequacy of the model [44–49]. Minimizing the weighted sum of square errors (5) corresponds to a Maximum Likelihood Estimate (MLE) of the parameters for the model in Eq (4).

The FIM for the model in Eq (4) is calculated as (6) The FIM is the inverse covariance matrix for the parameter estimates, so that the square root diagonal elements of the inverse FIM correspond to the one-standard deviation statistical uncertainties in the inferred parameters. It therefore follows that a well-conditioned FIM is necessary for accurate parameter inference.

Because mathematical models are always approximations, a model’s discrepancy from physical reality can always be improved by including additional physical mechanisms. This increased realism usually comes at the cost of increased complexity, often in the form of additional parameters, additional physical degrees of freedom, and computational cost. Therefore, a modeler is often faced with choosing from a hierarchy of models of increasing realism and growing complexity. Overly complex models increase the possibility of over-fitting data or over-explaining behavior while excessively simple models may lack all the relevant mechanisms. Effective models strike a balance between these two extremes.

To formalize this concept and facilitate later discussion, we introduce the concept of a sloppy system. We define a sloppy system as a physical system and a set of experimental protocols that can be approximated by a hierarchy of mechanistic, mathematical models of growing complexity that become sloppy in the limit of microscopic accuracy.

To make this definition more concrete, consider a model of a biological system. Of necessity, the model does not include every mechanism that is present in the true physical system. It is therefore possible to augment the model with additional mechanistic details, resulting in a more realistic, but more complex model. By repeating this process, one can generate a sequence of models of increasing realism that are better approximations to the actual physical system. The limit of this sequence converges on a model indistinguishable from the physical system, i.e., the limit of microscopic accuracy.

Naturally, this sequence of models will introduce many parameters that are unidentifiable. However, if the sequence of models is not just unidentifiable, but is also sloppy (i.e., as in the first column of Fig 1 rather than the third column), then we say that the system is sloppy.

As an aside, we have introduced the concepts of “sloppy systems” and “limits of microscopic accuracy” as useful abstractions. In practice, constructing more detailed mechanistic models may have a number of challenges. For example, it may be necessary to include parameters describing the experimental apparatus but that are separate from the biological system of interest. These potential complications are beyond the scope of this paper which instead uses these idealized concepts to explain our reported results.

The significance of a sloppy system is that models will always include marginal parameters so that there is always a trade off between model identifiability and model predictivity. For a fixed set of experimental protocols, there will always be some mechanistic details that are not entirely identifiable but still facilitate predictivity in the model.

The concept of a sloppy system brings together a number of concepts that are each well-known to modelers of complex systems. Considering multiple models of a single physical system is as old as science, but has recently been employed with new sophistication in the context of ensemble modeling [50] and multi-scale modeling [51–53]. The literature for parameter identifiability and the related concept of sloppiness was discussed in the introduction. By bringing these concepts together, we seek to explain the results of our EGFR simulations above and argue that sloppy systems pose unique challenges for predictive modeling.

We now demonstrate that the EGFR system is a sloppy system. Fig 6 shows the FIM eigenvalues for several models of EGFR signaling under the experimental conditions of Brown et al. [6]. Including more microscopic realism in the model requires additional parameters that render the model less identifiable. Note that this sequence of models was constructed in reverse–beginning from a sloppy model, irrelevant parameters were removed one at a time using the manifold boundary approximation method [17, 18]. Here we reinterpret this result as a demonstration of the existence of sloppy systems and emphasize the trade-off between mechanistic accuracy and model simplicity. In general, all the models that belong to a sloppy system cannot be ordered in a simple sequence as in Fig 6. There will often be a complex, hierarchical relationship among models similar to that described by the adjacency graphs in reference [54].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. An example of a sloppy system.

Observations of an EGFR signaling network can be explained by a model that is identifiable and not sloppy. The 18 parameter model has FIM eigenvalues that span fewer than 4 orders of magnitude and are all larger than one. By including additional mechanisms in the model (more parameters) the models become increasingly sloppy and less identifiable. The FIM eigenvalues ultimately span more than 16 orders of magnitude, leading to the large parameter uncertainties reported in reference [7].

https://doi.org/10.1371/journal.pcbi.1005227.g006

In order to account for errors due to ignoring marginal parameters, we need to refine the assumptions underlying Eqs (4)–(6). If the stochastic term in Eq (4) is dominated by experimental noise, then σ_i can be estimated by repeated observations of data point i, in which case σ_i scales like where n_i is the number of repeated observations of data i. In this case, σ_i becomes the standard deviation of the observations and is a measure of experimental reproducibility. Henceforth, we assume that σ_i denotes the experimental uncertainty and is known from experimental observations. We modify eq (4) to also include an error term δ_i due to approximations in the model (7) The model error term δ_i will be a type of “hyper-model” that accounts for and quantifies errors in the model in a phenomenological way without including additional mechanisms.

We adopt a simple hyper-model of the systematic error given by (8) where f is a hyper-parameter that will be estimated from the data, and is another Gaussian random variable with zero mean and standard deviation of one. We illustrate this concept geometrically in Fig 7. When this ansatz breaks down, it is an indication that relevant mechanisms are missing from the model, i.e., that the unfit data has structure that could be modeled and predicted.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Quantifying model error.

As in Fig 2, the model of interest forms a statistical manifold in data space, represented by the black dashed line. Another more realistic model also forms a statistical manifold of higher dimension (red surface). Experimental observations (blue dot) are generated by adding Gaussian noise of size σ to a “true” model (red dot). The least squares estimate is the point on the approximate model (black dot) nearest to the experimental observations. However, the distance from the best fit to the observed data has contributions from both the experimental noise and the model error.

https://doi.org/10.1371/journal.pcbi.1005227.g007

Care must be taken in the interpretation of . We have modeled the systematic error as a random variable. Unlike experimental noise, the size of this uncertainty cannot be decreased by repeated observations. Rather, the stochastic element in the model error represents the (unknown) approximations in the model. The relevant statistical ensemble is the set of all possible model refinements that could be made to correct the model shortcomings.

We have assumed that the model errors are uncorrelated among data points. We also assume that the model is likely to give worse predictions for data points that also have large experimental variation. These choices are convenient and constitute what is likely the simplest possible such hyper-model. More sophisticated models could be used, and the meta-problem of modeling the error in the model has been addressed in the context of uncertainty quantification [44–49]. In the present context, these assumptions will give us a simple way of estimating the error of the model from data. These assumptions will be valid provided that δ_i is small compared to experimental noise. We will now use our hyper-model to provide a criterion for including additional mechanisms in the model.

The negative log-likelihood for the model of Eqs (7) and (8) is (9) The best fit values for the parameters θ are unchanged by Eqs (7) and (8), and an unbiased estimate of f is given by (10) where χ² is the sum of squared error defined in Eq (5), M is the number of data points, and N are the number of parameters in the vector θ. As a practical matter, the model can be fit using Eq (5) as though there were no approximations in the model. This is due to our convenient choice of δ_i in Eq (8).

The augmented FIM is given by (11) where I in the first entry is the N × N FIM in Eq (6). The zeros in the off-digonal terms are 1 × N and N × 1 zero vectors. It follows that the parameter covariance matrix must be modified according to (12) where I is the FIM in Eq (6) and χ² is the best fit cost that minimizes Eq (5).

It is worth noting that the parameter f contains the same information as the likelihood function as is made explicit by Eq (10). Indeed, Eq (12) is a standard statistical formula for estimating parameter uncertainties in ordinary least squares regression in which the scale of the noise is unknown.

The standard deviation of the estimate in f is given by (13) We now consider how large an f can be acceptable. We seek a model’s whose approximations do not limit is predictive power. That is to say, the model error should be small compared to the experimental noise (f < 1) and unidentifiable from experimental observations. If the MLE of model error is f, then the statistical uncertainty in that estimate should satisfy δf ∼ f < 1. If this is the case, then the systematic error will be relatively small and not significantly limit the model’s predictive ability. This criterion gives (14) as an acceptable value for f.

EGFR Model Revisited

With this background, we can now revisit the EGFR model above. The experimental conditions of Apgar et al., included 7000 data points. If the approximate model were a good approximation, we would expect the fit to have a sum of squares error of approximately 7000 ± 84 where the range is one standard deviation of the Chi-squared distribution. However, fitting the artificial data typically led to a best fit error greater than 100,000 and was never less than 96,000. This error corresponds to an estimated value of f = 3.7 with δf = 0.03.

The statistical uncertainty in the model parameters is larger than expected by a factor for f = 3.7. The optimal experiments were designed to give less than 10% error in parameter estimates, so that modified uncertainty would now be less than 40%. Considering several parameters were initially unknown to a factor of a million, 40% appears to be a significant reduction. Although the parameter estimates remain somewhat constrained, the predictive power of the model is completely lost. This is because the effective error bars on the data are also larger by a factor of 3.8. In our simulations we assumed that the fractional activity levels of all proteins were measured to 10% accuracy. With 40% effective error bars, one standard deviation on either side of the mean covers almost the entire range of possible predictions.

In addition to having a large value for f, the ansatz of Eq (8) breaks down for the fitting the EGFR model. This is clearly seen by inspecting Fig 5. We speculate that it may be possible to rescue some of the predictive power of the model by implementing a more sophisticated hyper-model, such as introducing a separate f parameter for each time series or including phenomenological parameters to account for correlations in systematic errors. However, this possibility is beyond the scope of this work, but has been explored in the uncertainty quantification literature [44, 45, 47].

Mechanistic Mass-Action vs. Approximate Michaelis-Menten Kinetics

It is perhaps surprising that the approximate Michaelis-Menten model is inadequate even if Eq (3) is satisfied. However, one should remember that this condition was derived for a single enzyme-substrate reaction in isolation. One possible explanation of our results is that approximate Michaelis-Menten kinetics are not valid in a network. This explanation is problematic however, because the approximate Michaelis-Menten model had been used previously to fit real experimental data and make falsifiable predictions for new experiments. Indeed, in spite of its dubious status, the Michaelis-Menten approximation is often used with much success in many systems biology models. Therefore, while it is true that the Michaelis-Menten approximation is not generally valid, there is considerable evidence that it may sometimes be an safe approximation.

To complicate matters, we have forced Eq (3) to be satisfied by requiring K_M to be very large. Naively, one would expect this restriction to lead to K_m being unidentifiable. While this would also be true for measurements of a single reaction, it does not generalize to the network case, as our results demonstrate.

Furthermore, as the DNA repair results show, our results are not specific to the question of approximate Michaelis-Menten kinetics. Rather, we have shown that the general question of which physical details are necessary to include in a sloppy model can depend strongly, and in unexpected ways, on which combinations of experiments the model is to explain.

Implications for Optimal Experimental Design

A common use for optimal experimental design is model falsification. Demonstrating the shortcomings of model is hopefully accompanied by new insights into the system’s behavior. Since none of models we have considered here were considered “correct” in the reductionist sense, demonstrating that they are incomplete is not profound in itself. We suggested above that errors in the fit could be used to motivate new hypotheses about microscopic mechanisms. This possibility is beyond the scope of the current work that focuses on the implications for parameter estimation in sloppy models.

A potential alternative to experimental design for parameter estimation, is experimental design to constrain model predictions [28, 32]. Rather than constrain parameter estimates, one seeks to identify a small number of experimental observations that are controlled by the same few parameter combinations as the prediction one would like to make. In this approach, the model parameters remain sloppy, but the model may be predictive in spite of uncertainty about microscopic details.

It may be surprising that a model may be more predictive in the unidentifiable regime than in the identifiable regime. The predictivity of the an unidentifiable model is enabled by the narrow widths of the model manifold in Fig 2. The narrow widths guarantee that even infinite fluctuations in parameters do not correspond to large fluctuations in predictions. It has been suggested elsewhere that “sloppiness” can explain why models that make many uncontrolled approximations may be usefully predictive [9, 20]. Our results lend some support to this hypothesis; for our test cases removing sloppiness was always accompanied by a decrease in the predictive ability of the model.

Parameter Estimation in Sloppy Models

In previous work, sloppiness has been viewed as a challenge to be overcome or as a disease to be cured [5, 13, 15]. From this perspective the major challenge of sloppy models has been assumed to be the small eigenvalues of the FIM corresponding to practically unidentifiable parameter combinations. This in turn has led (incorrectly) to the conflation of sloppiness and practical unidentifiability. As we have argued here, the near uniform spacing of the eigenvalues (on a log scale) also pose unique challenges for parameter estimation because there is no clear cutoff between relevant and irrelevant mechanisms.

In order for an approximate model to be effective, it is important that the microscopic details ommited from the model be irrelevant, i.e., unidentifiable. When modeling systems for which all the relevant mechanisms are known, the validity of the model can usually be justified by small parameters, e.g., separated scales or distance from a critical point. The small parameter guarantees that the FIM eigenvalues for the irrelevant mechanisms are well-separated from those of the relevant mechanisms (e.g., column 3 in Fig 1). Some amount of unidentifiability in the physical system is therefore important for effective modeling.

For many complex systems, no such (known) small parameter exists and sloppy model analysis reveals that there is no sharp distinction between the relevant and irrelevant mechanisms. We speculate that in many cases the system (not just the model) is intrinsically sloppy because there is no intrinsic scale separation to suppress irrelevant mechanisms in the system. Therefore, a sequence of mechanistically more realistic models would have an eigenvalue structure closer to that in column 1 in Fig 1 rather than Fig 3. If that is the case, then one should not expect there to exist a mathematical model that can both be accurately calibrated and accurately predict the system behavior. There will always be several parameters that are marginal, i.e., not tightly constrained by data but are nevertheless necessary to explain the system behavior. In this case there is a fundamental limit to the efficacy of optimal experimental design: attempting to constrain the marginal parameters of a model of a sloppy system reduces the accuracy of the model and limits its predictive ability as we have seen.

Rather than posing a problem for parameter identification in models of complex systems, we argue here that sloppiness is important for successful modeling. Sloppy model analysis reveals that in many cases a behavior of interest is controlled by only a small number of parameter combinations. This observation has been used to explain why relatively simple models can make useful predictions. Indeed, it has been argued elsewhere that sloppiness may help explain why the world in its microscopic complexity is comprehensible at different scales [20]. Our results give credence to this position since removing sloppiness from a model reduced its predictive ability.

Another approach is to remove the sloppy parameters from the model. In principle, another simple model may exist whose parameters correspond to the few relevant parameter combinations in the sloppy model. Parameter estimation in such a model would be relatively straightforward. Recent advances in model reduction suggest that systematic construction of simple models from complex representations may be generally possible [17, 18, 20].

Relevant and Irrelevant Parameters in Sloppy Systems

In some branches of physics, the distinction between relevant, irrelevant, and marginal parameters is defined rigorously in terms of the stability of the collective behavior to microscopic variations in mechanistic details as measured by a renormalization group flow. In that context, relevant parameters correspond to degrees of freedom that must be tuned to achieve a behavior. In this work, we have used relevant and irrelevant less precisely as synonyms for identifiable and unidentifiable as measured by the FIM eigenvalues. This equivalence is reasonable because the the identifiable parameters are those that must be tuned to reproduce a behavior. The equivalence of these definitions was demonstrated in reference [12]. However, one of the hallmark features of sloppy models is the roughly uniform spacing of FIM eigenvalues, making it difficult to make a clear delineation between relevant and irrelevant parameters.

Lacking a clear cutoff between important and unimportant parameter directions means that some physical mechanisms may be either relevant or irrelevant depending on the experimental conditions. We have shown this explicitly for the two cases considered here. The model of Brown et al. [6] contained all of the relevant mechanisms (and many irrelevant ones) for explaining several experimental observations of an EGFR pathway. In contrast, it did not contain all of the relevant mechanisms for experimental conditions proposed by Apgar et al. [13]. Similarly, the LPL model is sufficient for modeling single radiation doses, while the RS model is necessary for modeling sequences of varied radiation doses and neither contains all the relevant mechanisms for modeling the experiments described in this work.

These results demonstrate the need for a theory of modeling and approximation that identifies which physical mechanisms are relevant for explaining different collective system behaviors. We have described two approaches that could be the beginnings of such a theory. First, we have introduced the concept of a sloppy system in which multiple models of varying complexity describe the same observations. Second, we have used a hyper-model to quantify the limitations of a model. Although, most of these ideas have existed in some form in the literature, the unique contribution of this work is synthesizing the concepts to explain why sloppy models pose unique challenges for system identification and why these problems are not shared by unidentifiable models that are not sloppy.

Because simple models are not complete, they cannot make accurate predictions for all experimental conditions. Of course, it is possible to extend a model by including more details in order to extend its range of validity. In principle, a single, monolithic model could accurately predict the outcome of all possible experiments. This possibility underlies the concept of a sloppy system. Microscopically complete models effectively act as numerical experiments and are a precursor to a more complete theory. In advancing to a more complete understanding of a system, we believe it is useful to consider multiple models of varying complexity and try to understand their limitations. Simultaneously considering multiple representations creates a rich and insightful theory into the mechanisms driving behavior that allow for abstraction and generalization. We believe that accounting for the approximations and context of a model are essential for successful modeling.

Conclusion

In this paper we have proposed a “sloppy system hypothesis.” We speculate that the prevalance of sloppy models in complex biological systems (and other areas of science) is not due to a limitation in the measurement structure of the system, but reflects a property intrinsic to the system itself. Because many complex systems lack an intrinsic scale separation (i.e., “small parameter” as discussed in the introduction), there is no mechanism whereby irrelevant details are necessarily suppressed in the model. Consequently, corrections to the mathematical model are relevant at all scales so that an accurate model will necessarily have several “marginal parameters” as in Fig 2. This hypothesis suggests that there is a fundamental limitation of optimal experimental design in sloppy systems due to these marginal parameters; attempting to constrain the marginal parameters of a model of a sloppy system reduces the accuracy of the model and limits its predictive ability. We have demonstrated this phenomenon on two complex biological systems, EGFR signaling and DNA repair.

Mathematical modeling in the face of structural uncertainty is a problem of growing importance across science [44–49]. Because mathematical models by their very nature are not exact replicas of physical processes, it is essential that they include the physical details relevant to the behavior of interest. In some branches of science, most notably several areas in physics, the equations governing some phenomena are sufficiently well-understood that numerical simulations come very near to being surrogates for real experiments. When this is the case, accurate parameter estimates reduce uncertainty in the model’s predictions. However, in many areas of complexity science, particularly for systems with fewer symmetries and less homogeneity, which physical details are relevant for explaining a particular behavior remains the theoretical bottleneck.

Our results suggest that there is a need for better understanding of and accounting for the approximations in complex models. In particular, optimal experimental design methods should limit their search space to those experiments for which the model is an accurate approximation. In spite of the growing documentation of microscopic biological mechanisms, it is difficult to predict how the errors introduced by a given approximation (such as the steady state approximation) will propogate to the predictions for the system’s collective behavior. In other words, it is difficult to know a priori which mechanisms are relevant for a particular behavior. We believe that better quantification of uncertainty will enable improved methods of experimental design and the development of accurate models for predicting behavior in complex systems.

Ethics Statement

Animals were maintained in an Association for Assessment and Accreditation of Laboratory Animal Care approved facility, and in accordance with current regulations of the United States Department of Agriculture and Department of Health and Human Services. The experimental protocol was approved by, and in accordance with, institutional guidelines established by the Institutional Animal Care and Use Committee.

UT MD Anderson Cancer Center

ACUF #: 00001061 RN00

Date Approved: 1/24/2014

Expiration Date: 1/2/2017

EGFR Models and Simulation

We use the model of EGFR signaling due to Brown et al. [7] formulated in terms of approximate Michaelis-Menten kinetics. We also constructed a similar model using mechanistic mass action kinetics. This replaces each approximate Michaelis-Menten reaction with two mass-action steps: E + S ⇌ ES → E + P. We first model each chemical reaction as an enzyme and substrate reversibly binding into an enzyme-substrate complex, and then dissociating to yield the original enzyme and the product. This gives four nonlinear ordinary differential equations (ODEs) for each enzyme substrate ppreaction, including one each for the changes in concentration of the enzyme, the substrate, the enzyme-substrate complex, and the resulting product. In total, modeling the EGFR network using the same topology as the Brown model by means of mechanistic mass action kinetics requires 54 independent, nonlinear ODES with 70 parameters. These equations are given in the supplement along with an sbml implementation of the mechanistic model and are available on github [60].

All models were simulated using in-house C, FORTRAN and python routines, including methods to automatically calculate parameter sensitivities, included in the supporting information. We use the approximate Michaelis-Menten model to simulate the original seven laboratory experiments performed by Brown et al. using parameter values from reference [7]. We then add random noise to results of this simulation and treat the results as if they were actual laboratory results for the experiments. Finally, we fit this data to the mechanistic mass-action model using the geodesic Levenberg-Marquardt algorithm [21, 61]. In order to help avoid complications in which the fit is prematurely stuck at manifold boundaries (as described in [21]), all fits were done with regularizing terms to keep parameters from drifting to infinite values. We use regularizing terms that take the form w_i(log x_i/x_i0)² for each parameter x_i. Fits were repeated for many different values of x_i0 (i.e., the point at which the regularization was centered) and weights w_i. Weights were varied over four orders of magnitude (0.01 to 100), and we observe that our final fits were robust to these choices, suggesting that our results are not an artifact of having converged to a local optimum.

Using the mechanistic model with parameters from fitting the Brown experiments, we simulate the five experimental conditions proposed by Apgar et al. [13] and add random noise. We then fit the approximate model to the these data as before.

Radiation Experiments

Experimental methods are the same as in reference [55]. C3Hf/KamLaw mice were exposed to whole body irradiation using 300 kVp X-rays at a dose rate of 1.84 Gy/min, and the number of viable jejunal crypts was determined using the microcolony assay. 14 Gy total dose was split into unequal first and second fractions separated by 4 h. Data were analyzed using the LQ model, the lethal potentially lethal (LPL) model, and a repair-saturation (RS) model.