## Abstract

The Wnt signalling pathway controls cell proliferation and differentiation, and its deregulation is implicated in different diseases including cancer. Learning how to manipulate this pathway could substantially contribute to the development of therapies. We developed a mathematical model describing the initial sequence of events in the Wnt pathway, from ligand binding to β-catenin accumulation, and the effects of inhibitors, such as sFRPs (secreted Frizzled-related proteins) and Dkk (Dickkopf). Model parameters were retrieved from experimental data reported previously. The model was retrospectively validated by accurately predicting the effects of Wnt3a and sFRP1 on β-catenin levels in two independent published experiments (*R ^{2}* between 0.63 and 0.91). Prospective validation was obtained by testing the model's accuracy in predicting the effect of Dkk1 on Wnt-induced β-catenin accumulation (

*R*≈0.94). Model simulations under different combinations of sFRP1 and Dkk1 predicted a clear synergistic effect of these two inhibitors on β-catenin accumulation, which may point towards a new treatment avenue. Our model allows precise calculation of the effect of inhibitors applied alone or in combination, and provides a flexible framework for identifying potential targets for intervention in the Wnt signalling pathway.

^{2}- computational modelling
- Dickkopf (Dkk)
- ordinary differential equation (ODE)
- secreted Frizzled-related protein (sFRP)
- signal transduction
- Wnt pathway

## INTRODUCTION

The Wnt signalling pathway controls mechanisms that direct cell proliferation, polarity and fate determination during embryonic development and in tissue homoeostasis. Numerous studies have linked mutations in the Wnt pathway to various diseases including cancer [1,2]. Many genes in the Wnt pathway, initially discovered as transiently functioning in development, have turned out to be oncogenes and tumour suppressors [3,4]. The special significance of Wnt signalling in stem and progenitor cells implies its role in regulating CSCs (cancer stem cells), which are thought to be associated with relapse, metastasis and drug resistance [5,6].

Significant effort is currently being invested worldwide in the development of therapeutic agents that function by manipulating the Wnt pathway, with particular focus on agents that may provide novel approaches to limiting tumour growth [3,6,7]. For example, researchers are currently developing anti-cancer therapeutic applications for regulatory proteins such as sFRP (secreted Frizzled-related protein) [8,9]. Since the Wnt pathway is involved primarily in embryogenesis and in tissue repair in adults, compounds that affect this pathway are not expected to cause significant toxicity [6].

There are at least 19 Wnt ligands and ten Frizzled receptors, activating at least three intracellular signalling pathways [4,10,11]. This vast network is influenced by a wide range of regulators including two main classes of extracellular inhibitors: those directly interacting with Wnt proteins [e.g. sFRP, WIF (WNT inhibitory factor 1)] and those binding to Wnt receptors or co-receptors [e.g. Dkk (Dickkopf) binding to LRP (low-density lipoprotein-receptor-related protein)] [12,13]. The complexity of interactions among ligands, antagonists and receptors makes the Wnt signalling pathway amenable to therapeutic intervention at many target points. It is not intuitively clear which component of the system is the best target for therapeutic intervention, or how such interventions should be designed in order to achieve the best clinical outcomes.

Mathematical models based on biological information are effective in improving the understanding of complex signalling pathways and their roles in disease control [14]. Such models are simplified to include only the main components of the signalling pathways, yet ensuring that the fundamental properties of the system are retrieved [15]. Theoretical and numerical analysis of such models may be used to predict system behaviour under various scenarios, and these predictions can then be compared with experimental data [15] {see [16] for an example of an application of this approach to the MAPK (mitogen-activated protein kinase) pathway}. However, as quantitative information on signal transduction pathways is rarely available, most mathematical models of signal transduction pathways are yet to be experimentally substantiated.

Among the Wnt-activated pathways, the canonical Wnt pathway is the best characterized. It regulates the transcriptional co-activator β-catenin, which controls the expression of specific target genes [4]. To study the dynamical behaviour of the intracellular components of this pathway, Lee et al. [17] formulated and parameterized a mathematical model using data from experiments in *Xenopus* extracts. Other researchers further analysed this model and extended it by adding pathway-related reactions (e.g. [18–21]; see review in [22]). These models focus on the intracellular steps of the canonical Wnt pathway. However, modelling the initial steps of this pathway, i.e. the binding of extracellular ligands and inhibitors to membrane receptors, is essential for elucidating the roles of different Wnt signal inhibitors. Such knowledge could pave the way to the development of innovative therapeutic interventions. For example, a relatively simple model for the Wnt and Notch signalling pathways suggested that exogenous Dkk1 is a potential modulator of stem cell fate decision, a prediction that was verified experimentally in mammary CSCs [23]. However, this model's simplicity, while rendering it analytically tractable, was an impediment to quantitatively evaluating the chemical reactions in these pathways [24].

In the present study, we developed a detailed mechanistic model for the extracellular and intracellular parts of the canonical Wnt pathway. The predictive ability of this general model was retrospectively validated using quantitative data from independent experiments, testing the effects of Wnt3a and sFRP on β-catenin accumulation. We also conducted experiments measuring inhibition by Dkk1, which prospectively validated model predictions of this effect. Simulations of the combined effects of sFRP1 and Dkk1 predicted synergism between these inhibitors.

## MATERIALS AND METHODS

### Mathematical model of the Wnt pathway

#### Basic assumptions and kinetic model

In our model, schematically described in Figure 1, the interactions among the Wnt ligand, its inhibitors and cell surface receptors are represented in detail, whereas the intracellular part is represented more generally. The model captures the following processes: (i) the pathway is activated by binding of Wnt ligand to the Frizzled receptor [25,26]; (ii) the resulting receptor–ligand complex may recruit an unoccupied LRP receptor and create a ternary complex consisting of Wnt, Frizzled and LRP [27–29]; and (iii) the latter complex transduces the signal inside the cell and interferes with the destruction cycle of β-catenin [30,31].

The intracellular level of β-catenin is regulated by a specific destruction complex comprising Axin, APC (adenomatous polyposis coli) and GSK3β (glycogen synthase kinase 3β). This complex binds β-catenin and causes its phosphorylation. Phosphorylated β-catenin dissociates from the destruction complex and is rapidly degraded [4].

Three assumptions of our model are consistent with those of Lee et al. [17]: (i) upon dissociation from phosphorylated β-catenin, the destruction complex may bind another β-catenin molecule; (ii) β-catenin is produced at a constant rate; and (iii) there exists an additional slow degradation path of β-catenin, independent of the destruction complex. Lee et al. [17] also assumed that the reactions assembling the destruction complex are at rapid equilibrium, with the exception of GSK3β binding. To further reduce the modelling complexity of the intracellular part, we assume that the destruction complex is at rapid equilibrium with all of its components. This implies that its total concentration is constant. Furthermore, we represent β-catenin phosphorylation, dissociation from the complex and degradation as a one-step process.

The regulation of β-catenin destruction by Wnt signalling is carried out via the intracellular domain of the LRP receptor. When bound in ternary complex with Wnt and Frizzled, the LRP receptor binds intracellular Axin [32,33]. For simplicity, our model assumes that the whole destruction complex is bound to the intracellular part of the activated LRP, in line with the above assumption of equilibrium between the complex and its components. We also assume that the binding of the destruction complex to the intracellular domain of LRP is reversible, and upon dissociation the receptor complex decomposes into its components [34]. The bound destruction complex cannot participate in the β-catenin destruction cycle. Hence, high Wnt concentration leads to the formation of active ternary Wnt–Frizzled–LRP complexes, which reduce β-catenin destruction, allowing its accumulation.

The model includes down-regulation of the pathway by sFRP, which competes with Frizzled for Wnt binding [35,36], and Dkk, which binds to LRP and abolishes formation of the ternary complex [13,37]. The total concentrations of all system components, except for β-catenin, are assumed to be constant over the time of interest.

#### Model equations

We translated the above description into ODEs (ordinary differential equations), assuming well-mixing and the law of mass action for all the protein interactions. The model describes the dynamics of 13 state variables, representing concentrations of the modelled pathway components by a system of seven ODEs. The system is closed by six conservation equations for the total concentrations of Wnt (*W _{T}*), destruction complex (

*C*), sFRP (

_{T}*S*) and Dkk (

_{T}*D*) and for the numbers of Frizzled (

_{T}*F*) and LRP (

_{T}*L*) receptors: The modelled variables represent the following: extracellular free Wnt (

_{T}*W*); Frizzled receptors free (

*F*) and bound to Wnt (

*F*); free LRP receptors (

_{W}*L*); ternary receptor complex Frizzled–Wnt–LRP (

*L*); intracellular destruction complex free (

_{F}*C*), bound by the ternary receptor complex (

*C*) and bound to β-catenin (

_{L}*C*); intracellular free β-catenin (

_{B}*B*); extracellular sFRP free (

*S*) and bound to Wnt (

*S*); extracellular free Dkk (

_{W}*D*); and LRP receptors bound to Dkk (

*L*).

_{D}Reaction rates are given by the coefficients *k*_{±i}, where *i* is the reaction step index as shown in Figure 1, and the sign corresponds to the reaction direction. Association reaction rates, *k*_{i}, have units of M^{−1·}s^{−1}, whereas dissociation rates have units of s^{−1}. Note that variables in different cell compartments have different units: intra- and extra-cellular proteins are measured in molar (M), whereas free or bound receptors in cell membrane are measured in the number per cell, receptor/cell. To compute reactions between components measured in different units, we introduced partition coefficients, *K _{su-ex}* and

*K*, translating numbers in cell membrane to extracellular and intracellular concentrations, respectively. These coefficients were determined as follows: and , where

_{su-in}*N*is the number of cells,

_{Cells}*V*is the experimental volume (both depending on experimental conditions),

_{exp}*V*is the intracellular compartment volume and

_{cell}*Av*is Avogadro's number. For the estimation of

*K*we assumed that experimental volume is a well-mixed external compartment. The estimation of

_{su-ex}*K*is based on a similar assumption for intracellular cytoplasmatic volume.

_{su-in}#### Computer implementation

The model was solved numerically by an ODE solver in MATLAB. Whenever some of the parameters were adjusted, it was performed by repetitive applications of a local search algorithm (trust region) in MATLAB, each time beginning with a random initial guess for the values of adjusted parameters; other parameters were set at literature-based values (see below).

### Data acquisition from the literature

To fine-tune and subsequently validate the model we used published experimental data on the time course of β-catenin accumulation induced by different Wnt3a concentrations (Figures 4C and 4D in [38]), as well as accumulation of β-catenin after a fixed period of stimulation by Wnt3a and inhibition of this effect by sFRP at different concentrations (Figures 1C and 2 in [39] and our previous work [40] respectively). All of the experiments were conducted in L cells, which do not express cadherins; hence, the cytoplasmatic β-catenin is free and not bound to the membrane. This makes these data suitable for comparison with our model, which does not include the binding of β-catenin to cadherins.

For all of the published measurements we extracted average values. We also report S.E.M. estimations, whenever available in the relevant publication source.

### Parameter evaluation

The initial estimation of parameter values and ranges was based on various literature sources (see the Supplementary online data at http://www.BiochemJ.org/bj/444/bj4440115add.htm). Then the data extracted from [38] were used to fine-tune the model for the Wnt effect. This was done in two steps: (i) adjustment of alternative subsets of model parameters to fit part of the experimental data (denoted the ‘partial training set’); and (ii) selection of the best-predictive parameter set, using the whole experimental data set.

We formed a partial training set by selecting three of the ten Wnt3a doses experimentally tested in [38] (0, 12.5 and 400 ng/ml). All of the time-course measurements of β-catenin accumulation for these doses were included in the partial training set.

At the first step, we considered several alternative choices of a subset of adjusted parameters, focusing on those parameters for which our initial estimation based on the literature was likely to be less reliable. Each subset of model parameters was adjusted by fitting model predictions to the partial training set, whereas other parameters were set at their initially estimated values. These values are reported in the Supplementary online data and the details of parameter adjustment are reported in the Appendix.

Subsequently, in order to select the best-predictive parameter set among the sets obtained in the first step, we simulated the application of all of the ten Wnt3a doses experimentally tested in [38], using each of these parameter sets, and compared the resulting model predictions to the data. To select the best parameter set, we used several statistical tests for correlation between observed and predicted values (see the Appendix). For example, we calculated the coefficient of determination, *R*^{2}, defined as , where *SS _{err}* is the sum of squares of differences between model-predicted and observed measurements: , where

*E*are the experimental data points,

_{i,j}*P*are the corresponding model predictions and

_{i,j}*M*is a matrix of weights (see the Appendix).

*SS*is the total sum of squares of the data, proportional to the sample variance: , where is the weighted mean of the observed data. This is the most general definition of

_{tot}*R*

^{2}; it implies that

*R*

^{2}≤1, and results that are closer to 1 indicate a better fit. Its value is not bounded from below, and when the prediction error is larger than data variance, it becomes negative, indicating a significant disagreement between the model and the data. This index was used also for comparing model predictions to experimental data during model validation. The best-predictive parameter set (Table 1) was used for further model simulations.

### Experiments testing Dkk1 inhibition of Wnt-induced β-catenin accumulation

#### Proteins

Recombinant mouse Wnt3a and Dkk1 were purchased from R&D Systems.

#### Assaying Wnt3a signalling by measuring β-catenin accumulation

The experiments were conducted under the same conditions reported in [40]. Mouse fibroblasts (L cells) were plated on 24-well plates at 250000 cells per well and grown overnight. They were incubated for 30 min in a serum-free medium, and then with Dkk1 at 0–10 nM for another 30 min. Wnt3a (0.5 nM final concentration) was then added. At 2 h later, the cells were lysed and β-catenin levels were detected by immunoblotting with anti-β-catenin mouse monoclonal antibody (BD Transduction Laboratories), and measured in Bio-Rad Laboratories ChemiDoc XRS. The membranes were re-blotted with mouse anti-β-actin antibody (Sigma). Numerical data of band intensity were obtained using Quantity One X software. β-Catenin band intensities were standardized with respect to β-actin band intensities. The densitometric intensity of experimental β-catenin bands was corrected by subtracting the negative control. These data were normalized to a scale where the relative densitometric value of β-catenin accumulation, induced by 0.5 nM of Wnt3a, 2 h after treatment was set at 100%. The results were expressed as the mean outcomes±S.E.M of three experiments analysed by ANOVA followed by a post-hoc multiple comparison test (using SPSS 16.0 software).

### Model validation

#### Predicting the effect of Wnt3a and sFRP

From two independent experimental studies [39,40], we extracted data on the inhibitory effect of sFRP (see the Data acquisition from the literature section). Both works report accumulation of β-catenin under several different concentrations of Wnt3a after 3 h [39] or 2 h [40]. In order to compare the results of the two experiments with model predictions, we separately scaled each one of them to the results of [38], using linear units scaling. The proportion coefficient between the measurement units of the experiment in [40] and those of [38] was found by minimizing the sum of squares of differences between β-catenin measurements in the two experiments at *t*=2 h, for the same concentrations of Wnt3a. This coefficient is equal to the proportion between the corresponding values of parameter λ, which translates the simulated β-catenin level to the experimental measurement units (see Table 1 and the Appendix). Similar scaling was performed for the data from [39], using the average between measurements at *t*=2 h and *t*=4 h in [38] for each Wnt3a concentration, since no measurements at *t*=3 h were available there. The obtained proportion coefficients were used to compute the values of λ for simulating each of these experiments.

To simulate the experiments testing the effect of Wnt3a, we used the same initial conditions as in the previously described simulation, and the experimental volume and cell number reported for each experiment. The concentrations of Wnt3a were set to 0.5 nM and 2.5 nM as reported for the experiments [40] and [39] respectively. The model was simulated over 2 or 3 h, corresponding to the duration of the respective experiment, and total computed β-catenin accumulation was compared with the experimental results.

Furthermore, the model was used to predict the inhibition effects of sFRP1 and sFRP2. To this end, we simulated the sFRP inhibition experiments, with *S _{T}* values ranging from 0 to 16 and from 0 to 500 nM, covering the range of experimental sFRP concentrations in [40] and [39] respectively. The inhibition of β-catenin accumulation was computed as the ratio between Wnt-induced β-catenin accumulation with and without sFRP. For this computation, the base level of β-catenin, i.e. that obtained for

*W*=

_{T}*S*=0, was subtracted from both values in the ratio, as done with the experimental results.

_{T}An additional adjustment of the rate parameters for sFRP2–Wnt binding was performed using the results of the sFRP2 inhibition experiment, which had been conducted with chicken sFRP2 [39] (whereas the binding coefficients in the model were initially estimated for mouse sFRP2 [40]). We used the search algorithm to determine the best-fit values for *k*_{1} and *k*_{−1}, whereas all of the other model parameters remained fixed.

#### Model validation using experimental results of inhibition by Dkk1

We used the model to predict Dkk1 inhibition of Wnt-induced β-catenin accumulation. The initial values were defined as above, and experimental volume and cell number were set in accordance with the Experiments testing Dkk1 inhibition of Wnt-induced β-catenin accumulation section. The values of *D _{T}* ranged from 0 to 10 nM. The inhibition was computed as the ratio between β-catenin levels with and without Dkk1, after 2.5 h. The baseline β-catenin level was subtracted from both values.

### Predicting the combined effect of sFRP and Dkk

We simulated the combined effect of sFRP1 and Dkk1 on Wnt-induced β-catenin accumulation, using the experimental volume and cell numbers as in the experiments testing Dkk1 inhibition of Wnt-induced β-catenin accumulation section. We checked different Wnt3a concentrations, in the range 0.05–5 nM. For each Wnt3a dose, we simulated the effects of adding different combinations of sFRP1 and Dkk1 concentrations in the experimentally relevant ranges of concentrations (0–40 nM for Dkk1 and 0–300 nM for sFRP1).

## RESULTS

### Model calibration

The mathematical model for the canonical Wnt pathway was constructed on the basis of published biological knowledge (see Figure 1 and the Materials and methods section), and calibrated by experimental data [38].

We found that adjusting as few as four parameters to a subset of the data, comprising measurements for three Wnt3a concentrations (0, 12.5 and 400 ng/ml), was sufficient to obtain a good fit to data for all other Wnt3a concentrations (Figure 2) with unbiased residual distribution (Figure 3). The adjustment of these four parameters resulted in the final set of parameters of our model (Table 1). Predictions generated using this parameter set, when compared with the complete data set [38], had a mean absolute error of 0.025 and coefficient of determination *R*^{2}≈0.965 (for definition see the Materials and methods section). These results demonstrate that our model accurately predicts the time course of Wnt-induced β-catenin accumulation.

### Model validation by independent experiments

#### Validation by predicting the effect of Wnt3a

The model predicted the effect of Wnt3a in independent experiments (Figure 4), with *R*^{2} values of 0.626 for the experiments that had employed Western blotting [40], and 0.908 for those that used ELISA as an assay [39]. Hannoush [38], whose data we used to fine-tune the model, also used ELISA, which may explain the better accuracy of model prediction for the experimental results of Galli et al. [39].

#### Validation by predicting sFRP inhibitory effect

Our model accurately predicts the sFRP1 effect on Wnt-induced β-catenin accumulation found in [40] and in [39] (Figures 5A and 5B respectively). The *R*^{2} values for predicted compared with the observed results are 0.893 and 0.911 respectively.

Prediction accuracy for the effect of sFRP2 was lower, *R*^{2}≈0.507 (Figure 6A) and *R*^{2}≈−0.65 (Figure 6B, solid line) respectively. In the latter case, the model clearly underestimated the experimental effect, as reflected by the negative coefficient of determination. However, the binding coefficients used in these model simulations were estimated for mouse sFRP2 [40], whereas the experiment shown in Figure 6(B) was conducted with chicken sFRP2 [39]. In this experimental work, it was also reported that commercially purchased mouse sFRP2 shows less potency than purified chicken sFRP2, and inhibits Wnt3a at doses comparable with those of sFRP1 [39]. Indeed, we found that increasing the binding rate to *k*_{1}≈8.89×10^{5} M^{−1}·s^{−1} and decreasing the dissociation rate to *k*_{−}_{1}≈7.8×10^{−5}·s^{−1} (both contributing to a decrease in the dissociation constant) render the model predictions of the experimental observations accurate (*R*^{2}≈0.914; see Figure 6B, broken line).

#### Prospective validation by experiments evaluating the inhibitory effect of Dkk1

The model predicted that low Dkk1 concentrations would inhibit Wnt3a-induced β-catenin accumulation (IC_{50}≈2.4 nM; see Figure 7A). Subsequently, our experimental results in L cells validated the model predictions with *R*^{2}≈0.944 (Figure 7); Dkk1 inhibited β-catenin accumulation in a dose-dependent manner with an IC_{50} of 3.2 nM.

### Model predictions of the combined effect of sFRP and Dkk

Figures 8(A) and 8(B) show simulation results of combined application of sFRP and Dkk in different concentrations, under 0.5 and 5 nM Wnt3a respectively. The effect of each inhibitor in the combination increased with its concentration; Dkk1 inhibited β-catenin to a specific level with lower concentrations than sFRP1 (also due to a lower *K*_{d} value). At a higher Wnt3a concentration, higher levels of the inhibitors are required to achieve the same relative inhibition.

These results were used to examine whether the combined effect of sFRP1 and Dkk1 is additive or synergistic. For this purpose, we created isoboles [41], curves that represent the set of sFRP1 and Dkk1 dose pairs that yield a specified effect. Each isobole in Figures 8(C) and 8(D) represents different combinations of sFRP1 and Dkk1 that inhibit β-catenin accumulation to a fixed level, in the presence of a given Wnt3a concentration. The convex form of the curves suggests a synergistic effect of the two inhibitors in all of the Dkk1 and sFRP1 concentrations tested, i.e. the effect of the combination is greater than additive, since an additive effect would have resulted in linear curves. Hence, lower concentrations of Dkk1 and sFRP1 can be combined to achieve significant inhibition of β-catenin. We also simulated application of the same concentration combinations, varying each of the model parameter values up to ±50%, and observed similar synergistic behaviour under all of the Dkk1, sFRP1 and Wnt3a concentrations tested (results not shown). This result implies that the observed synergistic effect is robust, probably imposed by the pathway structure.

Table 2 shows different combinations of sFRP1 and Dkk1 that are expected to inhibit Wnt-induced β-catenin accumulation to the same level, in the presence of a fixed Wnt3a concentration. The synergistic effect of Dkk1 and sFRP1 is clearly seen in Table 2 by comparing predicted concentrations of combined sFRP1 and Dkk1 (columns b, c and d) with the concentrations needed to achieve the same effect by each inhibitor alone (columns a and e). The ratio between sFRP1 and Dkk1 doses that separately achieve the same effect is termed the potency ratio. Column c shows the maximally synergistic combination of sFRP1 and Dkk1, i.e. the point where the sum of sFRP1 and Dkk1 doses, scaled by the potency ratio, is minimal. The ratio between Dkk1 and sFRP1 in the maximally synergistic combination increases when higher doses of Wnt3a are applied. Interestingly, this ratio is approximately the same as the potency ratio (Table 2).

## DISCUSSION

We employed mathematical modelling to investigate the Wnt canonical pathway and to provide a new mechanistic framework for this pathway, which can serve as a platform for discovering new anticancer treatments. Whereas previous mathematical models of the Wnt pathway have focused on its intracellular components (see review in [22]), our model details the initial sequence of events occurring along the pathway, from the ligands and inhibitors binding the membranal and extracellular components to β-catenin regulation. This allows exploration of quantitative effects of different extracellular pathway inhibitors, including therapeutic agents.

Our model accurately predicted the experimental results of the effects of Wnt and sFRP, obtained in different laboratories, by different assays and across a wide range of concentrations. In addition, model predictions concerning the effect of Dkk1 were prospectively validated experimentally. Note that the binding rates for these inhibitors, which are most influential on model predictions, were already fixed from the literature in the initial stage of parameter estimation. This demonstrates that our model is able to quantitatively predict the effects of inhibitors of the Wnt pathway.

For the effect of sFRP2 on Wnt-induced β-catenin accumulation, the model was less accurate (Figure 6). This can be explained by imprecise estimation of the rates of sFRP2 binding to Wnt3a. In one of the experiments [39] chicken sFRP2 was used, whereas the reaction rates assumed in our model were taken from mouse sFRP2 experiments [40]. Evaluating the chicken sFRP2–Wnt3a reaction rates by adjusting the model to the experimental data of Galli et al. [39], resulted in lower estimation for the dissociation constant. This concurs with the observation that mouse sFRP2 is less potent than chicken sFRP2 [39]. This demonstrates the applicability of our model for determining unknown specific parameters of Wnt pathway inhibitors (including potential drugs), when experimental data of their effect are available.

The Dkk1 dose–response relationship in our experiment varied from that recently observed by Bourhis et al. [42], with a difference of two orders of magnitude in the IC_{50} value. The reason underlying this discrepancy may be that in the latter study, the constant of dissociation (*K*_{D}) between Dkk1 and LRP was significantly higher than previously reported *K*_{D} values [13,37]. This difference may result from different protein expression systems and purification procedures.

Despite the abundance of published research on the canonical Wnt pathway, there are still parts of this pathway whose underlying mechanism is not clear. For example, our modelling hypothesis that a Wnt–Frizzled complex is created and subsequently binds to LRP is consistent with previous experimental results, but to the best of our knowledge had not been explicitly stated before. Experimental validation of our model lends support to the basic biological assumptions laid at its foundation, including this hypothesis.

Our general mechanistic model characterizes the canonical Wnt signalling pathway for different biological systems, and it is not tailored to a specific experimental setup. Hence, the model can be used to forecast results of various experiments, by input of a few parameters that characterize the specific experiment (e.g. cell number). New biological information, such as the binding affinities between specific Wnt ligands, their receptors and different inhibitors may also be easily incorporated into the model. In addition, new assumptions about the Wnt canonical pathway can be implemented in the mathematical model. For example, some studies suggest that Dishevelled polymerization induces large aggregates (signalosomes) of Wnt receptors and destruction complex (reviewed in [43]). This aggregation was best demonstrated when Wnt signalling proteins were overexpressed [44,45], and it still remains to be verified with endogenous Wnt signalling components. Once the information becomes available, this effect can be further incorporated into our model and its implications can be examined.

Since our model explicitly describes the extracellular part of the Wnt pathway, it can predict effects of various drugs that act on Wnt ligand–receptor interactions or on the modelled intracellular components. Therefore, it provides a flexible framework for planning an intervention in the Wnt signalling pathway and can also be used to predict the combined effect of drugs acting on both intra- and extra-cellular pathway components. The newly predicted synergism between sFRP1 and Dkk1 on Wnt-induced β-catenin accumulation should be experimentally verified. If experimentally supported, it may lead to the development of a new treatment modality for cancer and other diseases, namely, a combination of Wnt pathway inhibitors which are active through complementary mechanisms, like Dkk and sFRP.

In conclusion, our validated mathematical model of the initial stages of the canonical Wnt pathway is a useful tool for predicting the effects of different pathway inhibitors, alone or in combination. Our model provides a general framework for analysis of experimental data and systemic approach to investigating the Wnt pathway. Currently, we examine its use for studying the effects of different cancer-related mutations in the Wnt pathway and predicting effects of different pathway inhibitors or drugs on these mutants. The application of our model as a tool for testing the effect of possible therapeutic interventions should be studied further.

## AUTHOR CONTRIBUTION

Yuri Kogan, Karin Halevi-Tobias and Zvia Agur planned the work; Yuri Kogan and Karin Halevi-Tobias constructed the model; Yuri Kogan, Karin Halevi-Tobias and Gili Hochman estimated the model parameters; Yuri Kogan and Gili Hochman performed the numerical analysis; Luc Leyns planned the experiments and edited the paper prior to submission; Anna Baczmanska and Luc Leyns performed the experiments; Zvia Agur supervised the project; Yuri Kogan, Karin Halevi-Tobias, Gili Hochman and Zvia Agur wrote the paper.

## FUNDING

This work was partly supported by the NEST (New and emerging science and technologies) FP6 program [grant number 012930] and by the Chai Foundation.

## APPENDIX

### Parameter adjustment

Parameter adjustment to data from [38] was initially performed for several choices of parameter subsets. The considered subsets included four to seven parameters, e.g. {*F*_{T}, *k*_{3}, *k*_{6}, λ},{*k*_{3}, *k*_{5}, *k*_{6}, λ}, {*F*_{T}, *k*_{3}, *k*_{−5}, *k*_{6}, λ}, {*F*_{T}, *C*_{T}, *k*_{3}, *k*_{−5}, *k*_{6}, *k*_{−6}, λ}.

Given the vector *v* of the values of adjusted parameters, the model was simulated as follows. We set *S _{T}*=0 and

*D*=0, since in [38] no inhibitors were added, and used the reported experiment volume and cell number. For each Wnt3a concentration,

_{T}*W*, the model was initiated with free Wnt concentration (and

_{j}*W*) equal to

_{T}*W*, and other variables at their unique steady-state values assumed in the absence of Wnt. The model was numerically integrated from time

_{j}*t*=0 to

*t*=

*t*

_{10}=35 h. To compare the simulations with the experimental results, we interpreted the latter in terms of the model variables. We assumed that the fluorescence measurements reported in the experiments were in direct proportion to the total concentration of β-catenin, with proportion coefficient λ. Model predictions for β-catenin accumulation at time

*t*were computed as

_{i}*P*(

_{i,j}*v*)=λ[

*B*(

*t*)+

_{i}*C*(

_{B}*t*)], for a given Wnt3a concentration

_{i}*W*, and for parameter values

_{j}*v*. Here,

*B*and

*C*are the model variables representing, respectively, the concentration of free β-catenin and the concentration of β-catenin bound by the destruction complex. The value of λ was included in all subsets of adjusted parameters.

_{B}For each choice of a subset of adjusted parameters, we determined the parameter values by minimizing the GF (goal function), computed by comparing simulation results to the data in the partial training set. The data sampled from [38] were represented as a matrix *E*, where *E _{i,j}* is the β-catenin fluorescence intensity at time

*t*under concentration

_{i}*W*of Wnt3a, 1≤

_{j}*i, j*≤10. We also constructed a corresponding matrix

*M*of weights,

*M*∈{0,1,2} being the number of values available from graphical sampling of measurement

_{i,j}*E*. The partial training set was defined as {

_{i,j}*E*|

_{i,j}*i*=1,…,10;

*j*∈{1,5,10}}, corresponding to the doses

*W*

_{0}=0,

*W*

_{5}=12.5 and

*W*

_{10}=400 ng/ml).

We defined GF as the sum of squares of differences between model-predicted and observed measurements, with weights taken from matrix *M*:
where *J*={1,5,10}, and *P _{i,j}* are the corresponding model predictions under the parameter set

*v*. The values of the adjusted parameters that minimize the GF were found using the search algorithm with initial guesses randomly chosen within the relevant parameter ranges (reported in Supplementary online data).

### Tests used for the model calibration

The statistical tests used for evaluating the correlation between model-predicted and observed experimental values were: (i) mean residual between the simulation results and the experimental data points; (ii) average of the square of absolute and relative error per data point; (iii) the absolute value of absolute and relative error per data point; (iv) linear regression for predicted compared with observed data points; (v) linear regression as in (iv), forcing intersection with *y* axis to 0; (vi) coefficient of determination, *R*^{2} (defined in the Materials and methods section).

For tests (i)–(iii) smaller scores indicate better prediction; in test (iv) better prediction should give the intersection closer to zero and the slope and correlation coefficient closer to 1; in test (v) both the slope and correlation coefficient should be closer to 1; in test (vi) *R*^{2} should be closer to 1. When calculating average relative error in (ii) and (iii), the experimental points with values less than 0.005 (which is <1% of the maximal experimental result) were excluded from the dataset to prevent bias.

**Abbreviations:**
APC, adenomatous polyposis coli;
CSC, cancer stem cell;
Dkk, Dickkopf;
GF, goal function;
GSK3β, glycogen synthase kinase 3β;
LRP, low-density lipoprotein-receptor-related protein;
ODE, ordinary differential equation;
sFRP, secreted Frizzled-related protein

- © The Authors Journal compilation © 2012 Biochemical Society