## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Parameter-free methods distinguish Wnt pathway models and guide design of experiments

Edited* by Melanie H. Cobb, University of Texas Southwestern Medical Center, Dallas, TX, and approved January 23, 2015 (received for review August 28, 2014)

## Significance

The canonical Wnt/β-catenin signaling pathway is important for essential cellular functions and is implicated in many diseases. We introduce a new mathematical model that focuses on β-catenin degradation and protein shuttling between cellular compartments. We compare our model to others and show that all fit to time-dependent experimental data. To evade this parameter problem, we use algebraic methods and characterize model features that are independent of the choice of parameter values. We find that multiple responses to Wnt are feasible under certain conditions for the new model, but not for the others; moreover, we provide dependencies between species (variables) that inform future experiments and model discrimination. We also highlight the wide applicability of these tools across problems in systems biology.

## Abstract

The canonical Wnt signaling pathway, mediated by β-catenin, is crucially involved in development, adult stem cell tissue maintenance, and a host of diseases including cancer. We analyze existing mathematical models of Wnt and compare them to a new Wnt signaling model that targets spatial localization; our aim is to distinguish between the models and distill biological insight from them. Using Bayesian methods we infer parameters for each model from mammalian Wnt signaling data and find that all models can fit this time course. We appeal to algebraic methods (concepts from chemical reaction network theory and matroid theory) to analyze the models without recourse to specific parameter values. These approaches provide insight into aspects of Wnt regulation: the new model, via control of shuttling and degradation parameters, permits multiple stable steady states corresponding to stem-like vs. committed cell states in the differentiation hierarchy. Our analysis also identifies groups of variables that should be measured to fully characterize and discriminate between competing models, and thus serves as a guide for performing minimal experiments for model comparison.

The Wnt signaling pathway plays a key role in essential cellular processes ranging from proliferation and cell specification during development to adult stem cell maintenance and wound repair (1). Dysfunction of Wnt signaling is implicated in many pathological conditions, including degenerative diseases and cancer (2⇓–4). Despite many molecular advances, the pathway dynamics are still not well understood. Theoretical investigations of the Wnt/β-catenin pathway serve as testbeds for working hypotheses (5⇓⇓⇓⇓⇓⇓–12).

We focus on models of canonical Wnt pathway processes with the aim of elucidating mechanisms, predicting function, and identifying key pathway components in adult tissues, such as colonic crypts. We compare four preexisting ordinary differential equation models (5⇓⇓–8) and find, using injectivity theory, that for any given conditions and parameter values, none of the models is capable of multiple cellular responses.

In many tissues Wnt plays a crucial role in cell fate specification (3). At the base of colonic crypts, cells exist in a stem-like, proliferative phenotype in the presence of Wnt. As these cells’ progeny move up the crypt axis they enter a Wnt-low environment and change fate (perhaps reversibly), becoming differentiated, specialized gut cells (13). In neuronal and endocrinal tissues, Wnt/β-catenin data suggest cell fate plasticity under different environmental conditions (14, 15). Here, we introduce a new model motivated by experimental findings not described in previous models (16⇓–18) to investigate bistable switching in the Wnt pathway. We find the new model to be capable of multiple cellular responses; furthermore, our parameter-free techniques identify that molecular shuttling (between cytoplasm and nucleus) and degradation together may serve as a possible mechanism for governing bistability in the pathway, corresponding to, for example, a committed cell state and a stem-like cell state.

Comparison of models (and mechanisms) requires data; the type of comparison performed depends on the data at hand. If data show bistability (two distinct response states), then we could rule out all models that preclude bistability; however the converse is not true (a graded response may be compatible with all models). Experimental studies in *Xenopus* extracts have been performed to validate a model of Wnt signaling (5), with further pathway elucidation in refs. 19, 20; however, the parameters identified in these studies may differ markedly from those involved in mammalian Wnt signaling (21, 22). With the aim of discriminating between models, we present the five Wnt models under a unifying framework, with standardized notation to facilitate comparison. We fit parameters to recently published mammalian β-catenin signaling time-course data using Bayesian inference (22) and find that all of the studied models can describe the data well, demonstrating that additional data are required to compare models.

To determine which sets of protein species should be measured for carrying out a comparison between models and data, we introduce matroid theory to systems biology. A matroid is a combinatorial structure from mathematics, and in our case, it provides all of the steady-state invariants (23, 24) that have minimal sets of variables. The algebraic matroid associated with the steady-state ideal determines specific sets of species that should be measured to perform model discrimination without knowledge of parameter values. We demonstrate this parameter-free analysis with experimental data for two Wnt models.

In the next section, we introduce the previous models and new shuttle model. We perform injectivity–multistability analysis and classify the shuttle model as the only one capable of multistability. Next we infer the parameters of five competing models for time-course β-catenin data, revealing that all of the models fit the data. Finally we introduce algebraic matroids to inform experimental design for discriminating between models and data.

## Models

Over the past decade, Lee et al.’s seminal model of canonical Wnt signaling (5) has spawned many variants. Briefly, the underlying biology of the pathway that these models describe is as follows (25): Wnt binds to cell-surface receptors that transduce a signal via a multistep process involving Dishevelled (Dsh) to the so-called destruction complex (DC), which contains forms of Axin, adenomatous polyposis coli (APC), and glycogen synthase kinase (3) (GSK-3). In the absence of a Wnt signal, the DC actively degrades β-catenin––which is being continually synthesized in the cell––by phosphorylating it and marking it for proteasomal degradation. Following Wnt stimulation, degradation of β-catenin is inhibited through phosphorylation of DC members, leading to accumulation in the cytoplasm of free β-catenin, which is able to translocate to the nucleus where it can form a complex with T-cell factor (TCF) and lymphoid-enhancing factor proteins and influence the transcription of target genes. These are cell-type specific, although genes controlling self-renewal and proliferation are commonly regulated across many cell types (1).

We include the core processes as well as the following in the proposed shuttle model (Fig. 1*A* and *SI Appendix*):

*i*) spatial localization of shuttling components Axin, APC, GSK-3, and Dsh, the importance of which has been reported for each species (17, 26⇓⇓⇓–30);*ii*) an alternative degradation mechanism whereby β-catenin is degraded while still bound to active DC and sequestered but not degraded by inactive DC (16);*iii*) catalysis of the reverse reaction by phosphatase (P) that converts DC from inactive to active form by dephosphorylating members of the DC (18, 31, 32).

The behavior of four other published models is analyzed and compared with that of the shuttle model (5⇓⇓–8). Fig. 1*B* summarizes the distinguishing qualitative features of each model; full model descriptions, using a standardized notation that differs from the authors’ originals, are summarized by a composite model in *SI Appendix*.

## Results and Discussion

Wnt signaling interaction networks are polynomial systems whose steady-state solutions are defined by sets of algebraic equations for the species’ concentrations; this opens up avenues for parameter-free analysis, as we show here.

### Parameter-Free Analysis of Wnt Models I: Multistability.

We are interested in determining whether or not a given model can produce multiple positive stable responses (states). Standard approaches from dynamical systems [e.g., bifurcation and singularity theory (33⇓–35)] are useful for small systems or when we have knowledge of the parameters, however, for systems of more than a few free variables (the shuttle model has 19 species and 31 parameters), such approaches become infeasible. To overcome this, we apply theory developed for chemical reactions to Wnt pathway models; this is particularly helpful for determining whether multiple states are possible without the need for parameter values or sampling. There are various conventions in chemical reaction network theory (CRNT) for describing the number of positive steady states; we use the following terminology:

*i*) injectivity implies at most one steady state;*ii*) multistationarity is the capacity for multiple steady states;*iii*) multistability is the capacity for multiple stable (accessible) steady states.

We test the injectivity of each model following graph-theoretic or Jacobian-based approaches used in CRNT (36⇓⇓⇓–40). We find that only the Schmitz et al. (7) and the shuttle model fail injectivity and exhibit multistationarity. Further analysis reveals that the Schmitz et al. model is capable of at most two steady states, only one of which is stable (see *SI Appendix* for proofs). Whereas all of the previous models possess at most one positive stable steady state for any choice of the parameter values and conserved species’ concentrations, the shuttle model has the capacity for multistationarity and multistability. We find that when three species shuttle (e.g., Dsh, inactive DC, and β-catenin), the model exhibits two stable states; we proceed to analyze this version of the shuttle model. Previous mathematical studies have proven shuttling across compartments is a mechanism for multistability (41, 42).

For each model, we have a minimal collection *SI Appendix*). Negating these gives necessary, but not sufficient, conditions on the parameters for multistationarity, which depend explicitly on the degradation and shuttling rate constants (Fig. 2*B*). Because the parameter values are unknown, we compute an illustrative bifurcation diagram for a specific parameter set satisfying these necessary conditions. These diagrams are similar to dose–response curves in an experimental setting, and they demonstrate different behavior as the shuttling and degradation rate constants are varied (Fig. 2 *A* and *C*).

Within a certain parameter region, we can either observe a graded response or switching and hysteresis behavior (Fig. 2*C*). The hysteresis loop shown by the black (stable state), blue (switch to committed state), and red (threshold switch to stem-like state) arrows enables switching at different thresholds between two stable steady states. Whereas over short timescales bi- and monostable behavior is indistinguishable (Fig. 2*D*), at steady state these differences emerge. In the bistable regime the low level of gene transcription is associated with a committed cell state and the high level with a stem-like phenotype over long time periods (Fig. 2*D*). As the value of a parameter, for example β-catenin shuttling into the cytoplasm *A*).

If qualitative data showed a clear bistable switch, then the shuttle model would be the best model. However, given quantitative rather than qualitative data, how can we compare models?

### Wnt Model Comparison via Parameter Inference.

Where competing models describe the same biological processes, one can perform parameter inference or model selection; such methods have been applied to a variety of problems in systems biology, ranging from cancer modeling to population genetics (43, 44).

Inferring parameters from data via Bayesian analysis provides the posterior probability distribution over the parameters, from which more information can be gleaned than by point estimates alone. In Fig. 3*B*, we demonstrate the Bayesian inference procedure by considering a 2D subset of the parameter space (rates of β-catenin synthesis and β-catenin degradation). This panel shows how over successive iterations we can home in on the most probable region of parameter space given the data.

Each of the Wnt models has a different number of parameters. In an attempt to compare the models fairly and to reduce the size of the parameter space that we are searching, for each model we choose to fix all of the parameters (at estimated or arbitrary values) except for three. These three are allowed to vary and are used to fit the model to the Wnt/β-catenin time-course data recently published in ref. 22; see *Materials and Methods* for details. We chose the free parameters based on their point of influence on the pathway, targeting parameters with direct or near-direct influence on β-catenin dynamics (Fig. 3*A*); the fits we obtained after performing inference are shown in Fig. 3*C*. We see that even with only three degrees of freedom, good parameter fits are obtained for all of the models. Studying the posterior for each model reveals relationships between parameters: high β-catenin production and low β-catenin degradation rates are favored across models, but β-catenin-TCF binding rates vary considerably between models.

The disparity between model complexity and data availability prevents us from choosing between models based on model selection analysis. The problem could be addressed by simplifying models or collecting more data [additionally, experimental design influences Bayesian model selection results (45)]. Here, we proceed to use parameter-free methods to help guide experiments for model discrimination.

### Parameter-Free Analysis of Wnt Models II: Matroids.

Instead of classifying the feasible behaviors of the whole system, we can use the finer structure of a model to derive relations in each part of the system, for example, the concentrations of species in a chemical reaction network at steady state. The matroid of a model is a list of

*i*) the subsets of species that are related, and*ii*) the subsets of species whose concentrations are unrelated.

A matroid is a set with a notion of independence for its subsets. The classic example of a matroid is an arrangement of vectors. Suppose *A*. A set of vectors is called “dependent” if it is linearly dependent, i.e., there is a set of scalars, not all zero, such that multiplying by the vectors and adding them together results in the zero vector. Here, *A*). The matroid does not remember the coordinates of the vectors, only whether any subset is dependent or independent. The matroid construction in our discussion uses algebraic independence instead of linear independence. Explicitly, if there is a polynomial relationship that a collection of species satisfies at steady state, they are considered a dependent set. Note that this only considers relationships at steady state, where the possible species concentrations describe an algebraic variety. The independent and dependent sets of molecular players in each model may help compare models, guide experiments, and possibly reject models as described throughout this section.

We calculate the matroid of five Wnt signaling pathway models (four shown in Fig. 4*B*). Each model has a rank *r*, which dictates the number of species from the full set whose concentrations can be independently specified; taking measurements of *r* independent species determines (in terms of parameters) the values for all other species. Circuits are minimal dependent sets of species––they become important when we consider model discrimination. A matroid can be represented pictorially by point arrangements: the set of species labeling a point has rank 1, the set on a line has rank 2, the set on a two-dimensional plane has rank 3, and so on (Fig. 4*C*). Any two species labeling the same point are algebraically related, as are any three species on a line, any four species in a plane, etc.

We describe how the matroid of a model is computed in more detail in *SI Appendix*; the input is the polynomial ideal of steady-state relations and the output is a list of all circuits with their polynomial relations. Strictly speaking, a circuit is defined as a set of variables. However, in this application, we record the polynomial relations, because these are the support-minimal steady-state invariants. One approach involves computing a Gröbner basis for every elimination ordering, a feasible although lengthy computation for small systems. An alternative uses linear algebra to pinpoint the sets of variables appearing in invariants; then, it uses this information in conjunction with elimination or numerical algebraic geometry software to find polynomials (see ref. 46 for more detail). This approach is only now being implemented because algebraic matroids have only recently been adopted for applications, e.g., low-rank matrix completion (47).

The results of the matroid calculations (Fig. 4*B*) prompt biological insight; for illustration, we analyze the van Leeuwen et al. model (12). In this model, five species (called “loops”) can be determined from just the parameters. Among the others, any pair not including

Unlike the other models, the solution set for the shuttle model has two irreducible components (loosely, proper subsets that should be considered separately). The matroids for the two irreducible components both have rank 5. The number of components for a given model is determined by its algebraic structure; the number of components and real positive steady states are not analogous––multiple steady states may appear even when we have only one component. As described above, if we want to know all species concentrations in the shuttle model, five measurements must be made and these measurements should be chosen to be independent. For example, measuring TCF (*T*) and any four species not lying in the same plane in Fig. 4*D* would determine all species concentrations.

The minimal sets of species from the matroid can also be used to study part of the system. The partial information (relations on a subset of species) obtained from the matroid can be used as a self-consistency check between data and competing models and, in this way, serves as a method to rule out models. Model discrimination based on steady-state invariants has been performed for one specific ordering, so by including the matroid we can recover all possible steady-state invariants with different variables (full set of circuit polynomials, as defined in ref. 48), including conserved quantities. For example, Tan et al. (22) measured average β-catenin in the cytoplasm (*X*) and the nucleus *X* and *X* and *δ*). From the data in ref. 22, we can test its compatibility with the model via a parameter-free method as described in ref. 24. Briefly, the method tests whether the data satisfy the Schmitz et al. circuit polynomial by checking whether there exist

Model compatibility is determined by computing the coplanarity error *X*. The null hypothesis that the model is compatible with the data can be rejected when the coplanarity error (normalized smallest singular value) is greater than a statistical bound as described in ref. 24 and *SI Appendix*, which is determined by the Gaussian measurement noise in the data and the invariant structure.

Before we test the models with data from ref. 22, we simulate data from both the shuttle and Schmitz et al. models. We draw random parameters from a lognormal distribution and then simulate 100 replicate measurements of *X* and *SI Appendix*). By consulting the matroid of the shuttle model, we find that *X* and

## Conclusions

There is a wealth of mathematical and experimental research on Wnt signaling, aimed at understanding the pathway well enough to target Wnt-implicated diseases. There are two significant challenges to overcome. The disparity between models and data that we have highlighted via Bayesian inference prevents us from constraining parameter values in a manner that often helps to elucidate mechanisms and predict function. The second challenge is the gap between in vitro and in vivo studies, and the corresponding differences in parameter values. This is supported by evidence on the variation in parameter estimates between naive and crowded (physiological) in vitro molecular experiments (49, 50). To gain insight into these complex systems described by complicated models, we must evade this parameter problem.

Parameter-free approaches can provide additional information about the β-catenin/Wnt pathway. Based on injectivity–multistationarity analysis, we find that the shuttle model predicts the possibility of a regulatory switch, acting early in the cell fate determination pathway. Other systems have also reported early checkpoints in cell fate signaling in activation of apoptosis through receptor–ligand binding (51, 52). We identify the possibility of important roles for spatial localization and degradation in cell fate switching. In the Erk pathway, it is also seen that either localization (via shuttling) or degradation by apoptosis is crucial for bistable switching, both mathematically and experimentally (42, 53⇓–55). We report for the first time, to our knowledge, that a combination of these processes governs the dynamical regime.

By computing the algebraic matroid of different Wnt models, we can characterize the dependencies between species. The matroid results enable us to guide experiments (which species to measure) to discriminate between models with data, all of the while not requiring parameter values. The mechanisms for bistability identified above are of course not exclusive and one could imagine many other models exhibiting bistability via different mechanisms. This scenario could prompt previously unidentified insight using our model discrimination framework, now with multiple bistable models.

Given the current (and growing) complexity of models across a wide range of topics, tools such as those demonstrated here offer strong potential for testing models and for predictions to be made. In addition, we provide possible directions for future experimentation to narrow the gap between data and models, and, through our predictions, help to unravel the workings of the intricate and essential Wnt pathway.

## Materials and Methods

### Bayesian Inference.

Model selection in systems biology can be performed using Bayesian inference (56). Here we perform parameter inference for model selection using approximate Bayesian computation, which forgoes evaluation of the likelihood function and instead calculates the (here Euclidean) distance between model and data (57), implemented in the ABC Sys-Bio package (58). For each model we compare the total free β-catenin level (in some cases addition of two species) with the data provided by ref. 22.

### Multistability.

Determining whether a model is capable of multiple responses can be tested using injectivity. The injectivity of each model was determined using CRNT toolbox (59); for those that were not injective (Schmitz et al. and shuttle), we computed the determinant of the Jacobian following (42), and analyzed the sign of the coefficients in Mathematica (Version 9.0; Wolfram Research). Bifurcation diagrams were computed using Oscill8 (available at oscill8.sourceforge.net/doc) and visualized with MATLAB (R2013a; The MathWorks).

### Matroids.

Matroid computation is performed by structured variable elimination on eligible polynomial systems. This is carried out using symbolic algebra software Macaulay2 (60) with the aid of packages presented in ref. 46. Code is available at math.berkeley.edu/∼zhrosen/matroids.html. When the set of steady-state solutions has multiple irreducible components, the matroid was computed for each to assess the independence structure in each regime. Isolated points in the solution set were not analyzed, as the matroid is trivial. Model discrimination was performed in Sage (available at cloud.sagemath.com) following the method presented in ref. 24.

## Acknowledgments

We thank A. Burgess and C. Wee Tan for data and discussions about Wnt signaling. We also thank the anonymous reviewers as well as T. Dale, E. Feliu, A. Fletcher, K. Ho, P. K. Maini, E. O’Neill, A. Shiu, and B. Sturmfels for helpful discussions and/or comments on the manuscript. H.A.H. gratefully acknowledges funding from Engineering and Physical Sciences Research Council Fellowship EP/K041096/1 and the American Institute of Mathematics. Z.R. and H.A.H. acknowledge funding from Royal Society International Exchanges Scheme 2014/R1 IE140219. A.L.M. and H.M.B. acknowledge funding from the Human Frontiers in Science Program (RGP0039/2011). A.L.M., H.M.B., and H.A.H. acknowledge funding from King Abdullah University of Science and Technology KUK-C1-013-04.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: harrington{at}maths.ox.ac.uk.

Author contributions: H.M.B. and H.A.H. designed research; A.L.M., Z.R., and H.A.H. performed research; A.L.M., Z.R., and H.A.H. analyzed data; and A.L.M., Z.R., H.M.B., and H.A.H. wrote the paper.

The authors declare no conflict of interest.

↵*This Direct Submission article had a prearranged editor.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1416655112/-/DCSupplemental.

## References

- ↵
- ↵.
- Polakis P

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Li X,
- Yost HJ,
- Virshup DM,
- Seeling JM

- ↵
- ↵.
- Hernández AR,
- Klein AM,
- Kirschner MW

*β*-catenin specify the sites of Wnt control. Science 338(6112):1337–1340 - ↵
- ↵
- ↵
- ↵.
- Harrington HA,
- Ho KL,
- Thorne T,
- Stumpf MPH

- ↵
- ↵.
- Franca-Koh J,
- Yeo M,
- Fraser E,
- Young N,
- Dale TC

- ↵.
- Wiechens N,
- Heinle K,
- Englmeier L,
- Schohl A,
- Fagotto F

- ↵.
- Cong F,
- Varmus H

- ↵.
- Henderson BR,
- Fagotto F

- ↵
- ↵.
- Luo W, et al.

*β*-catenin degradation complex. EMBO J 26(6):1511–1521 - ↵.
- Strovel ET,
- Wu D,
- Sussman DJ

- ↵.
- Guckenheimer J,
- Holmes P

- ↵.
- Kuznetsov YA

- ↵.
- Golubitsky M,
- Schaeffer D

- ↵
- ↵
- ↵.
- Craciun G,
- Tang Y,
- Feinberg M

- ↵
- ↵
- ↵
- ↵
- ↵.
- MacLean AL,
- Filippi S,
- Stumpf MPH

- ↵.
- Jiang R,
- Tavaré S,
- Marjoram P

- ↵
- ↵.
- Rosen Z

- ↵.
- Király FJ,
- Theran L,
- Tomioka R

- ↵.
- Király FJ,
- Rosen Z,
- Theran L

- ↵.
- Aoki K,
- Yamada M,
- Kunida K,
- Yasuda S,
- Matsuda M

- ↵
- ↵
- ↵
- ↵.
- Li B,
- Dou QP

- ↵
- ↵.
- Michailovici I, et al.

- ↵
- ↵
- ↵.
- Liepe J, et al.

- ↵.
- Ellison P,
- Feinberg M,
- Ji H

- ↵.
- Grayson DR,
- Stillman ME

## Citation Manager Formats

### More Articles of This Classification

### Physical Sciences

### Applied Mathematics

### Biological Sciences

### Related Content

- No related articles found.