Thierry Castel

Christophe Lecomte

Yves Richard

Annabelle Larmure

Victor Enduro

Lucy Terese

Modelling the regional evolution of frost damage on a winter crop with warming in a temperate climate Simulation results in Burgundy-Franche-Comté

The modalities of climate warming over Europe and France allow to distinguish two temperature periods: before and after 1987/1988 (Brulebois et al. 2015 J. Hydrol.). It offers the opportunity to assess climate warming on agrosystems.

Perspectives Objectives & Methods

Winter Pea (Pisum sativum L)

The results show that the account of the winter crop traits and the modalities of warming jointly matter to better understand how winter climate warming affects changes in the frost risk evolution. This is of primary importance to define new winter crop ideotypes.

An automated mapping based on a kriging with external drift (KED) method has been used to interpolate maximum and minimum temperatures.

Regional frost stress evolution

Climate data (Burgundy-Franche-Comté -Northeast France) recorded by the Météo-France Station Network (MFSN) were daily interpolated on a 12-km grid.

Frost resistance (°C) Acclimation rate (days)

The interpolated climate data were used to feed a frost stress crop model validated for pea (Lecomte et al. 2003 Agronomie;Castel et al. 2017 OCL).

Sowing dates: october 1 st , 10 th , 20 th and 30 th In average, the number and the intensity of frost stress events evolve in opposition between both periods (1961-1987 and 1988-2015): the number of events increases, while their intensity decreases, especially for variety with a weak frost resistance and a prolonged acclimation rate. The resulting evolution of frost stress (number of events x intensity) is a decrease (Castel et al. 2017 OCL).

Average frost stress evolution

  Sowing date: october 30 th P1.114


Konstantinos Koutroumpas
email: konstantinos.koutroumpas@ecp.fr

Paolo Ballarini

Irene Votsi

Paul-Henry Cournède

Bayesian parameter estimation for the Wnt pathway: An infinite mixture models approach

come

Bayesian parameter estimation for the Wnt pathway: an infinite mixture models approach

Introduction

Mathematical modeling has become an integral part of modern biological research. With the onset of Systems Biology, mathematical models have been widely used to formally describe complex biological systems. Models enable the use of mathematical analysis and computer simulations not only to understand but also to predict and recently even control the behavior of biochemical processes [START_REF] Milias-Argeitis | In silico feedback for in vivo regulation of a gene expression circuit[END_REF][START_REF] Uhlendorf | Long-term model predictive control of gene expression at the population and single-cell levels[END_REF]. A challenging task during model building is the estimation of its parameters from experimental data. The equations describing the dynamics of a system depend on several parameters. For instance, chemical reactions are described by rate parameters, which capture how quickly a protein is synthesized or degraded or how fast a protein moves from one location to another in the cell. While some parameters can be measured using biochemical studies, chemical-kinetic models also include a substantial number of rate parameters that must be inferred from experimental data.

Parameter estimation approaches can be classified into two categories: frequentist and Bayesian. In the frequentist approach a point, or interval, estimation of model parameters is provided by maximizing the likelihood function L(θ θ θ) = p(D o |θ θ θ), naturally defined as the probability that the observed data D o have been generated by parameters θ θ θ.

Bayesian methods, on the other hand, have been proven a powerful alternative to frequentist methods, one in which parameters are treated as random variables with distributions attached to them (i.e. posterior distributions p(θ θ θ|D o )). The overall goal in this case is to estimate the parameters' posterior distribution based on: 1) prior beliefs/knowledge on the unknown parameters, in form of the prior distribution p(θ θ θ), and 2) knowledge coming from the experimental data, in form of the likelihood function L(θ). By relying on probability distributions, Bayesian methods do not require large sample approximations to provide one with interval estimates but rather allow one to use prior beliefs on the unknown parameters, as expressed through the prior distribution p(θ θ θ). In other words, in a Bayesian context, the prior beliefs are updated by the observed data D o via the likelihood function p(D o |θ θ θ), and the posterior distribution is determined up to a normalized constant as follows p(θ θ θ|D o ) ∝ p(θ θ θ)p(D o |θ θ θ).

In several cases the calculation of the likelihood function L(θ θ θ) is difficult or even impossible, precluding the use of conventional likelihood-based inferential techniques. However, simulating from the likelihood is generally straightforward, even if obtaining a reliable numerical/functional representation of the model is not possible. In such settings, likelihood-free computation or approximate Bayesian computation (ABC) [START_REF] Tavaré | Inferring coalescence times from DNA sequence data[END_REF][START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF][START_REF] Marjoram | Markov chain Monte Carlo without likelihoods[END_REF][START_REF] Sisson | Sequential monte carlo without likelihoods[END_REF] can be employed. ABC methods circumvent the explicit evaluation of the likelihood using a simulation-based approximation. The original ABC algorithm jointly simulates θ θ θ ∼ p(θ θ θ) and D s ∼ p(D s |θ θ θ ) and accept θ θ θ if and only if the simulated dataset D s is close to the observed dataset.

The simple ABC scheme suffers from the same shortcamings as other rejection samplers: most of the samples are drawn from regions of parameter space, which cannot give rise to simulation outputs that resemple data. Therefore a number of computational schemes have been proposed to enhance the efficiency of ABC algorithms. ABC coupled to SMC samplers [START_REF] Toni | Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems[END_REF] aim to sample sequentially from a sequence of distributions, which increasingly resemble the target posterior. The successive distributions are constructed by sampling parameters from the previous samples, and perturbing them according to a transition kernel [START_REF] Del Moral | Sequential monte carlo samplers[END_REF]. The choice of the transition kernel has a significant impact on how accurate the approximate posterior is [START_REF] Filippi | On optimality of kernels for approximate bayesian computation using sequential monte carlo[END_REF]. In most cases a random-walk kernel is used to locally perturb a sample [START_REF] Givens | Local adaptive importance sampling for multivariate densities with strong nonlinear relationships[END_REF][START_REF] Filippi | On optimality of kernels for approximate bayesian computation using sequential monte carlo[END_REF]. However, the choice of the transition kernel in an ABC context remains an open question, since the kernel influence is not formally understood and can be problem-dependent.

Here we revisit the problem of optimal transition kernel design using infinite mixture models. Finite mixture models have been previously used as importance densities in Adaptive Importance Sampling approaches [START_REF] Cappé | Adaptive importance sampling in general mixture classes[END_REF][START_REF] Wraith | Estimation of cosmological parameters using adaptive importance sampling[END_REF][START_REF] Cornuet | Adaptive multiple importance sampling[END_REF], with however a difficulty in choosing the proper number of elementary distributions in the mixture. We thus propose the use of infinite mixture models as successive important sampling distributions in the SMC framework, which will prove to bring very appealing results as detailed below. We formally describe how infinite mixture models and more precisely Dirichlet Process Mixtures (DPM) can be integrated in the ABC-SMC algorithm. We show that the ABC-SMC with DPM kernels is significantly more efficient in approximating complex, multimodal distributions and we illustrate its use in a model of the Wnt signalling pathway using real experimental data. Given the wide use of the ABC-SMC algorithm in parameter estimation and model selection and the improved performance of the proposed variant, we expect that it will significantly impact future modeling efforts in Systems biology.

The remainder of the paper is organized as follows. A short description of the ABC and ABC-SMC algorithms is provided in Section 2, followed by a description of the Dirichlet process Gaussian mixture model. We then present how DPGMM can be used as a transition kernel in the ABC-SMC framework. In Section 3 some biological background on the Wnt signalling pathway, the biological process we are interested in, is given accompanied by a mathematical model describing the main chemical reactions in the pathway. The methodology is applied both to a toy example and the Wnt pathway model in Section 4. Finally, some concluding remarks and future extensions are discussed in Section 5.

Methods

Approximate Bayesian Computation

Approximate Bayesian Computation (ABC) methods [START_REF] Tavaré | Inferring coalescence times from DNA sequence data[END_REF][START_REF] Pritchard | Population growth of human Y chromosomes: a study of Y chromosome microsatellites[END_REF][START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF] replace the calculation of the likelihood with direct sampling from the posterior p(θ θ θ|D) ∝ p(D|θ θ θ)p(θ θ θ) through simulation. Given the prior distribution p(θ θ θ) on the parameter space, a candidate parameter vector θ θ θ * is sampled from p(θ θ θ) and a simulated data set D s is generated from the likelihood function p(D|θ θ θ * ). The parameter vector is accepted if the simulated data set is close to the observed data set D o according to a distance function d (d(D o , D s ) ≤ ). The distance function detrermines the discrepancy between the two datasets and is a tolerance value. If the data are too intricate or complicated it is common to replace a comparison of the observed and simulated data by a comparison of suitable summary statistics. An ABC rejection algorithm was used in [START_REF] Pritchard | Population growth of human Y chromosomes: a study of Y chromosome microsatellites[END_REF] for a model of the variation in human Y chromosome and its basic scheme follows Algorithm 1.

Input: Prior distribution p(θ θ θ), observed data set D o , likelihood function p(D|θ θ θ), number of accepted particles N , tolerance

for i = 1 → N do Sample θ θ θ * from p(θ θ θ). Simulate a data set D s ∼ p(D|θ θ θ). while d(D o , D s ) > do Sample θ θ θ * from p(θ θ θ). Simulate a data set D s ∼ p(D|θ θ θ). end Set θ θ θ i = θ θ θ * end Algorithm 1: ABC rejection algorithm
The output of the ABC rejection algorithm is a sample of parameter vectors from the distribution p(θ θ θ|d(D o , D s ) ≤ ). In the case that is sufficiently small, the distribution p(θ θ θ|d(D o , D s ) ≤ ) is considered a good approximation of the posterior distribution p(θ θ θ|D o ). The acceptance rate of the ABC rejection algorithm depends on the proposal distribution p(θ θ θ) and it can become very low for a diffuse prior. [START_REF] Marjoram | Markov chain Monte Carlo without likelihoods[END_REF] proposed an ABC method based on Markov Chain Monte Carlo (ABC-MCMC) to overcome the low acceptance rate problem. Although the algorithm is guaranteed to converge to the approximate posterior distribution p(θ θ θ|d(D o , D s ) ≤ ) while substantially improving the acceptance rate over the ABC rejection algorithm [START_REF] Marjoram | Markov chain Monte Carlo without likelihoods[END_REF], the samples in the generated Markov Chain are highly correlated and the chain may stuck in areas of low probability for long periods of time [START_REF] Sisson | Sequential monte carlo without likelihoods[END_REF].

Sequential Monte Carlo for ABC

Algorithms combining ABC and Sequential Monte Carlo sampling have been proposed [START_REF] Sisson | Sequential monte carlo without likelihoods[END_REF][START_REF] Toni | Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems[END_REF] to avoid some of the disadvantages of the ABC-MCMC methods. To improve acceptance rates due to the mismatch between the proposal and the target distribution, sequential sampling methods [START_REF] Del Moral | Sequential monte carlo samplers[END_REF] define a sequence of intermediate distributions p 1 , . . . , p T that evolve gradually from a proposal distribution to the target distribution. The intermediate distributions are approximated by a set of random samples, the so-called particles. Starting with a distribution p 1 that is easy to approximate using the prior p, a number of particles, {θ 1 θ 1 θ 1 , . . . , θ N θ N θ N }, are sampled from the prior and targeting p 1 . The particles are then moved using a Markov kernel. It corresponds to sampling from an importance sampling distribution and the new particles are thus weighted with importance sampling weights [START_REF] Toni | Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems[END_REF], taking into account how well a particle conforms to the next target distribution and the importance sampling distribution used to generate the samples. In the ABC framework the intermediate distributions are constructed by considering different tolerance values, i.e. p i = p(θ θ θ|d(D o , D s ) ≤ i ) where 1 > 2 > . . . > T ≥ 0. In this study an ABC-SMC algorithm [START_REF] Toni | Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems[END_REF] based on sequential importance sampling [START_REF] Del Moral | Sequential monte carlo samplers[END_REF]) is used and its scheme for a deterministic system is shown in Algorithm 2.

Apart from the increased efficiency in comparison to the ABC rejection algorithm, ABC-SMC does not stuck in low probability regions and the particles are not correlated, which is the case with ABC-MCMC [START_REF] Sisson | Sequential monte carlo without likelihoods[END_REF]. Moreover, given that a population of particles is used, the method is more efficient when complex posteriors are considered [START_REF] Sisson | Sequential monte carlo without likelihoods[END_REF]. However, the performance of the algorithm clearly depends on the selected tolerance sequence. Fast decreasing tolerances can result in low acceptance rates and poor performance, while if they decrease slowly the algorithm becomes computationally intensive [START_REF] Del Moral | An adaptive sequential Monte Carlo method for approximate Bayesian computation[END_REF]. For the optimal selection of the tolerance sequence, Del [START_REF] Del Moral | An adaptive sequential Monte Carlo method for approximate Bayesian computation[END_REF] proposed an adaptive ABC-SMC algorithm. The kernels K j that are used to draw new particles starting from the previous ones also affect ABC-SMC efficiency. Being easy to implement, local random-walk moves (Gaussian or uniform) are usually employed [START_REF] Sisson | Sequential monte carlo without likelihoods[END_REF][START_REF] Toni | Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems[END_REF]. The parameters of the kernels are determined by some statistics on the previous sample population. Such approaches can be problematic especially in the case of complex multimodal posteriors [START_REF] Filippi | On optimality of kernels for approximate bayesian computation using sequential monte carlo[END_REF]. [START_REF] Filippi | On optimality of kernels for approximate bayesian computation using sequential monte carlo[END_REF] have derived optimality criteria based on the Kullback-Leibler divergence between a target distribution and its estimate. They have also showed that for complicated posterior distributions, locally adapted kernels that are designed taking into account the K-nearest neighbors for each particle tend to show the best performance. The main drawback of this approach is that the number K has to be fixed a Input: Prior distribution p(θ θ θ), an observed data set D o , a likelihood function p(D|θ θ θ), number of accepted particles N , a tolerance schedule 1 > . . . > T ≥ 0

for j = 1 → T do if j == 1 then for i = 1 → N do Sample θ θ θ * from p(θ θ θ). Simulate a data set D s ∼ p(D|θ * θ * θ * ). while d(D o , D s ) > j do Sample θ θ θ * from p(θ θ θ). Simulate a data set D s ∼ p(D|θ * θ * θ * ). end Set θ θ θ j i = θ θ θ * ; w j i = 1 end Normalize the weights w i j else for i = 1 → N do Sample θ * θ * θ * from N i=1 w j-1 i K(θ θ θ|θ θ θ j-1 i ).; Simulate a data set D s ∼ p(D|θ * θ * θ * ). while d(D o , D s ) > j do Sample θ * θ * θ * from N i=1 w j-1 i K(θ θ θ|θ θ θ j-1 i ).; Simulate a data set D s ∼ p(D|θ * θ * θ * ). end Set θ θ θ j i = θ θ θ * ; w j i ∝ p(θ θ θ j i ) N i=1 w j-1 i K(θ θ θ j i |θ θ θ j-1 i )
.

end Normalize the weights w j i . end end Algorithm 2: ABC-SMC algorithm priori and its value affects the exploration of the parameter space [START_REF] Filippi | On optimality of kernels for approximate bayesian computation using sequential monte carlo[END_REF]. For this reason, we propose to use a Dirichlet Process Gaussian Mixture Model to generate the proposal distributions at each iteration of the SMC-ABC scheme.

Dirichlet Process Gaussian Mixture Model

Gaussian mixtures (GM) are a class of parametric mixture models that are used for density estimation [START_REF] Mclachlan | Finite mixture models[END_REF]. [START_REF] Cappé | Adaptive importance sampling in general mixture classes[END_REF], [START_REF] Wraith | Estimation of cosmological parameters using adaptive importance sampling[END_REF] and [START_REF] Cornuet | Adaptive multiple importance sampling[END_REF] have proposed the use of GM's as importance densities in Adaptive Importance Sampling approaches. As with the K-nearest neighbors approach in [START_REF] Filippi | On optimality of kernels for approximate bayesian computation using sequential monte carlo[END_REF], the number of components has to be preselected and in most approaches model selection is used to select the optimal number of components. A workaround to this problem is the use of a non-parametric Infinite Mixture model that assumes that the data come from a mixture of an infinite number of distributions. In the following we consider the Dirichlet Process Mixture (DPM) model [START_REF] Antoniak | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF] and more precisely the Dirichlet Process Gaussian Mixture model (DPGMM) which has been extensively used for density estimation and data clustering [START_REF] Müller | Bayesian curve fitting using multivariate normal mixtures[END_REF][START_REF] Rasmussen | The Infinite Gaussian Mixture Model[END_REF][START_REF] Escobar | Estimating normal means with a Dirichlet process prior[END_REF][START_REF] Maceachern | Estimating normal means with a conjugate style Dirichlet process prior[END_REF][START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF]. Starting with a finite GM, we can obtain a DPGMM by taking the limit as the number of components goes to infinity.

A GM with K components can be written as:

p(x x x| (π 1 , µ µ µ 1 , Σ 1 ) , . . . , (π K , µ µ µ K , Σ K )) = K c=1 π c N (µ µ µ c , Σ c ) (1)
where π c , µ µ µ c and Σ c are respectively the mixing proportion, the mean vector and the covariance matrix for component c respectively. If we assume a Dirichlet prior on the mixing proportions and a prior distribution H on the mean and the covariance, the model can be written:

x x x i |z i , (µ 1 , Σ 1 ) , . . . (µ K , Σ K ) ∼ N (µ µ µ zi , Σ zi ) z i |π 1 , . . . , π K ∼ Cat(π 1 , . . . , π K ) π 1 , . . . , π K |α ∼ Dir(α/K, . . . , α/K) (µ µ µ i , Σ i ) ∼ H (2)
where z i is the indicator variable and corresponds to the component to which the observation x i belongs. Given a set of observations {x 1 , x 2 , . . . , x N } the probability of an observation to belong to one of the K components is [START_REF] Görür | Dirichlet process gaussian mixture models: Choice of the base distribution[END_REF]:

p(z l = j|z z z \l , α) = C \l,j + α/K N -1 + α (3)
where C \l,k = i=1,...,N i =l δ k (z i ) is the number of observations except l that belong to component j. Taking the limit as K goes to infinity the probability that an observation belongs to a non-empty component (i.e components with C \l,j > 0) is:

p(z l = j|z z z \l , α) = C \l,j N -1 + α (4)
and the probability to belong to one of the infinitely many empty components is:

p(z l = z k ∀ k = l|z z z \l , α) = α N -1 + α . (5)
Probabilities (4-5) are the same with the probabilities in a Chinese restaurant process that is used to define a Dirichlet process. Hence, the infinite limit of a GMM (2) is a DPGMM [START_REF] Görür | Dirichlet process gaussian mixture models: Choice of the base distribution[END_REF] which can be compactly written:

x x x i | (µ µ µ i , Σ i ) ∼ N (µ µ µ i , Σ i ) (µ µ µ i , Σ i ) ∼ G G ∼ DP(H, α). ( 6
)
where DP(H, α) is the Dirichlet process with concentration α and base distribution H.

The base distribution H is the mean of the DP, which means that the DP draws distributions "around" H as a normal distribution draws real numbers around its mean. A strong property of DP is that the distributions drawn are almost surely discrete even if the base distribution is continuous. The concentration parameter α specifies how strong the discretization is. For an α close to 0, the distributions are concentrated on a single value, while as α goes to ∞ the realizations are close to H. The component parameters are then sampled from the distribution drawn from the DP. Both the choice of the base distribution and the value of the concentration parameter are important for the performance of the model. While the concentration parameter can sometimes be inferred from the data, it is harder to decide on the base distribution. In most cases the base distribution is chosen so as to facilitate the inference of the DPGMM model.

DPGMM estimation

DPM models are usually inferred using a MCMC algorithm [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF]. Most methods are based on Gibbs sampling, where variables (parameters, hyperparameters and indicator variables) are updated in turn by sampling from their posteriors conditional on the other variables. In the case that the base distribution is the conjugate prior of the probability distribution of the observations generated from a component, Gaussian in our case, a collapsed Gibbs sampler [START_REF] Maceachern | Estimating normal means with a conjugate style Dirichlet process prior[END_REF][START_REF] Sudderth | Graphical models for visual object recognition and tracking[END_REF] is more efficient than explicitly sampling parameters. In this case, the state of the Markov chain constructed from this sampler consists of the indicator variables only that are updated by iteratively reassigning each observation to a new component based on the assignments of the remaining observations. Consider a set of observations x x x = {x 1 , x 2 , . . . , x N } from a GM (1) with unknown number of components K. The probability of an observation l to belong to a component j is:

p(z l = j|z z z \l , x x x, α, κ, θ, ν, ∆) = p(z l = j|z z z \l , α) (7) p(x l |z l = j, z z z \l , x x x \l , κ, θ, ν, ∆)
where p(z l = j|z z z \l , α) is given by ( 4) or (5) and p(x l |z l = j, z z z \l , x x x \l , κ, θ, ν, ∆) can be seen as the likelihood of x l to have been generated from the same Gaussian with the observations that already belong to the component j. The latter can be approximated by a Gaussian with mean and covariance computed by the observations belonging to the component j [START_REF] Sudderth | Graphical models for visual object recognition and tracking[END_REF]. The Gibbs sampler starts by assigning all observations to a single component, and iteratively reassigns each observation to a new component based on the assignments of the remaining observations according to (7).

In the case of DPGMM a natural selection for H is the Normal-Inverse-Wishart distribution, which is the conjugate prior of normal distribution. The covariance matrix Σ is drawn from an inverse-Wishart distribution:

W -1 (Σ|ν, ∆) ∝ |Σ| (-ν+d+1 2 ) exp - 1 2 tr(ν∆Σ -1 ) (8)
where ∆ is the covariance parameter, ν the degrees of freedom and d the dimension. Conditioned on Σ the mean µ µ µ is drawn from a Normal distribution

N (θ, Σ/κ). (9)
A detailed description of the collapsed Gibbs sampler that was used for this study can be found in Section 2.5.3, Algorithm 2.3 in [START_REF] Sudderth | Graphical models for visual object recognition and tracking[END_REF]. For the concentration parameter α an inverse gamma prior was used as proposed in [START_REF] Rasmussen | The Infinite Gaussian Mixture Model[END_REF].

ABC-SMC using DPGMM kernels

In the ABC-SMC framework, DPGMM can be used to fit a mixture on the samples of each population and new samples can be drawn from the identified Gaussian components. More precisely, assume that

Θ j-1 Θ j-1 Θ j-1 = {θ j-1 k θ j-1 k θ j-1 k } N k=1 are the samples accepted during iteration j -1 of the ABC-SMC algorithm, W j-1 W j-1 W j-1 = {w j-1 k } N k=1
the corresponding weights, and z j-1 z j-1 z j-1 = {z j-1 k } N k=1 the indicator variables estimated using the previously described Gibbs sampler. Samples for the next iteration (lines 16 and 19 in Algorithm 2) can be directly drawn from:

K j-1 c=1 πc N (µ µ µ j-1 c , Σ j-1 c ) ( 10
)
where K j-1 is the number of non-empty components of the DPGMM and πc = It is important to note that the sampling weights are also considered in the calculation of the parameters of each component, i.e. the mean and the covariance. The mean of a component c is given by the weighted mean:

µ j-1 c µ j-1 c µ j-1 c = k∈C j-1 c w j-1 k θ t-1 k θ t-1 k θ t-1 k (12)
and the covariance matrix by the weighted covariance matrix:

Σ j-1 c = V j-1 c V j-1 c 2 -V j-1 c ( 13
) k∈C j-1 c w j-1 k θ j-1 k θ j-1 k θ j-1 k -µ j-1 c µ j-1 c µ j-1 c T θ j-1 k θ j-1 k θ j-1 k -µ j-1 c µ j-1 c µ j-1 c where V j-1 c and V j-1 c
are the sum of the weights and the sum of the squares of the weights.

V j-1 c = k∈C j-1 c w j-1 k (14) V j-1 c = k∈C j-1 c w j-1 k 2 (15)
3 The Wnt pathway

Wnt pathways are a family of intracellular signalling transduction pathways known for playing a key role in cell's development (e.g. cell proliferation, stem cell maintenance, differentiation, cardiac development).

They are the subject of intensive studies both in relation to cancer development as well as to embryonic development pathologies [START_REF] Fleiss | Microglial wnt signaling inhibition promotes microglia activation and oligodendrocyte maturation blockade[END_REF]. Two families of Wnt pathways have been, so far, characterized according to their dependency on the multifunctional β-catenin protein [START_REF] Kühl | An updated overview on wnt signaling pathways: A prelude for more[END_REF].

Canonical Wnt pathways (Wnt/β-catenin) regulate gene transcription via β-catenin mediated transduction of Wnt signal, whereas non-canonical Wnt pathways, like Wnt/Ca 2+ [START_REF] De | Wnt/ca2+ signaling pathway: a brief overview[END_REF], are involved in the regulation of a cell's physical characteristics, i.e. the cell shape, and its functioning. All Wnt pathways share the same activation mechanism: the signal is triggered through the binding of a Wnt protein to a receptor (of the Frizzled family) on the cell's membrane, and propagated inside the cell through the cytoplasmic dishevelled (Dsh) proteins. In this paper we consider a model of the canonical Wnt pathway.

Wnt/β-catenin signalling (basics). Signal transduction along the Wnt/β-catenin pathway is initiated when extracellular Wnt proteins bind to receptors of the Frizzled family and to coreceptors such as the lipoprotein receptor-related protein (LRP5/6). It triggers a cascade of events that involve a number of actors both in the cytoplasm and in the nucleus of the cell. The core species known to be involved in the pathway include (Table 1): the Dishvelled (Dsh) protein, the β-catenin protein, the Axin protein, the adenomatous polyposis coli (APC), the glycogen synthase kinase 3 (GSK-3β), the multimeric protein complex APC/GSK3β/Axin (given by complexation of APC, Axin and GSK-3β) also called the destruction complex (DC), and at nucleic level, the T-Cell factor (TCF) and the lymphocyte enhancing factor (LEF) which yield gene expression.

In the absence of extracellular Wnt molecules, β-catenin is destructed by the degradation complex, which forms around the activated form of Axin (AxinP), in the cytosol. When Wnt binds to LRP at the membrane, the assembly of the destruction complex is inhibited, leading to an accumulation of β-catenin in the cytosol. Consecutively, β-catenin is relocated to the nucleus where it activates the transcription of genes including the gene encoding for Axin protein.

Models of canonincal Wnt pathway. Many computational models of the Wnt/β-catenin pathway have been proposed [START_REF] Mazemondet | Elucidating the Sources of β-Catenin Dynamics in Human Neural Progenitor Cells[END_REF][START_REF] Maclean | Mathematical and statistical techniques for systems medicine: The wnt signaling pathway as a case study[END_REF][START_REF] Jensen | A wnt oscillator model for somitogenesis[END_REF][START_REF] Schmitz | Analysing the impact of nucleo-cytoplasmic shuttling of β-catenin and its antagonists APC, Axin and GSK3 on Wnt/β-catenin signalling[END_REF][START_REF] Lee | The roles of apc and axin derived from experimental and theoretical analysis of the wnt pathway[END_REF][START_REF] Tan | Analysis of wnt signaling β-catenin spatial dynamics in hek293t cells[END_REF]. They may be distinguished with respect to the biological features they capture, for example whether they consider the compartmentalization of the cell, or the granularity with which molecular activities are represented [START_REF] Maclean | Mathematical and statistical techniques for systems medicine: The wnt signaling pathway as a case study[END_REF]. Some of the models [START_REF] Tan | Analysis of wnt signaling β-catenin spatial dynamics in hek293t cells[END_REF][START_REF] Schmitz | Analysing the impact of nucleo-cytoplasmic shuttling of β-catenin and its antagonists APC, Axin and GSK3 on Wnt/β-catenin signalling[END_REF][START_REF] Mazemondet | Elucidating the Sources of β-Catenin Dynamics in Human Neural Progenitor Cells[END_REF] account for the cytoplasm/nucleus relocation of some of the species. In [START_REF] Tan | Analysis of wnt signaling β-catenin spatial dynamics in hek293t cells[END_REF] the dynamics of β-catenin compartmental diffusion have been investigated through a double-compartment ODEs model of Wnt-pathway which has been validated through time series measurements on Humans Epithelial Kidney (HEK293T) cells. In [START_REF] Jensen | A wnt oscillator model for somitogenesis[END_REF] In the context of this paper, we also focus on the Lee model and present some numerical experiments (Section 4) performed on it. We proceed below with a detailed description of the Lee model.

The Lee model.

The model consists of 6 core species (Dsh, Axin, β-catenin, TCF, APC and GSK3) and their basic interactions, (i.e. synthesis, degradation, phosphorylation, reversible complexations). The overall dynamics of the model are described by 15 mass action kinetics ordinary differential equations (16-29). In the [START_REF] Lee | The roles of apc and axin derived from experimental and theoretical analysis of the wnt pathway[END_REF] model, the existence of intracellular compartments (i.e. cytosol, nucleus) is abstracted away, hence all species are assumed to be contained in a single cellular compartment. Furthermore, the Wnt signal is not explicitly represented, rather it is assimilated to the state of the dishevelled (Dsh) protein (active Dsh corresponding to the presence of Wnt signal on the membrane receptors, inactive Dsh corresponding to its absence).

In absence of Wnt signal (Dsh i ), the DC (Axin * /APC * /GSK3) is free to form through a multistep reaction, consisting of: first the reversible dimerization of Axin with APC (APC/Axin, binding at rate proportional to α 8 , unbinding at rate proportional to α 9 ), then the reversible formation of the Axin/APC/GSK3 trimer (binding at rate α 6 , unbinding at rate proportional to α 7 , if Dsh is inactive, or to α 5 if Dsh is active), followed by the reversible phosphorylation of Axin and APC parts of the Axin/APC/GSK3 trimer to obtain the actual DC (Axin * /APC * /GSK3). The DC can then reversibly sequester free β-catenin molecules (sequestration at rate proportional to α 10 , liberation at rate proportional to α 11 ). β-catenin sequestration by the DC eventually leads to β-catenin degradation (at rate proportional to α 14 ) through an intermediate, irreversible, phosphorylation step (phosporylation at rate proportional to α 12 ). Notice that the decomplexation of the unphospohorylated form of DC (i.e. the Axin/APC/GSK3 complex) is regulated by the active form of Dsh (i.e. Dsh a ): in presence of Dsh a the decomplexation of Axin/APC/GSK3 is enhanced by a rate which is proportional to the level of Wnt signal (i.e. the amount of Dsh a ) as well as to the constant rate α 5 . The model also accounts for constant production and degradation of both β-catenin (produced at constant rate α 15 , degraded at rate proportional to α 16 ) and Axin (produced at constant rate α 17 , degraded at rate proportional to α 18 ). Activation/inactivation of Dsh (corresponding to the binding/unbinding of Wnt molecules to the cell's membrane) are assumed to occur at rate proportional to α 1 , respectively α 2 .

R 1 : Dsh i α 1 -- α 2 Dsh a R 2 : APC * /Axin * /GSK3 α 4 -- α 3 APC/Axin/GSK3 R 3 : APC/Axin + GSK3 α 6 ------------------ α 5 [Dsha]+α 7 APC/Axin/GSK3 R 4 : APC + Axin α 8 -- α 9 APC/Axin R 5 : APC * /Axin * /GSK3 + β-cat α 12 ---- α 11 APC * /Axin * /GSK3/β-cat R 6 : APC * /Axin * /GSK3/β-cat α 12 APC * /Axin * /GSK3/β-cat * R 7 : APC * /Axin * /GSK3/β-cat * α 13 APC/Axin/GSK3/ + β-cat * R 8 : β-cat * α 14 ∅ R 9 : ∅ α 15 β-cat R 10 : β-cat α 16 ∅ R 11 : ∅ α 17 Axin R 12 : Axin α 18 ∅ R 13 : β-cat + T CF α 19 ---- α 20 β-cat/T CF

A compartmentalized Lee model

An important feature of the [START_REF] Lee | The roles of apc and axin derived from experimental and theoretical analysis of the wnt pathway[END_REF] model to study the role of the Wnt pathway in diseases is that it focuses on the assembly of the DC from its constituent parts. While available experimental data are not adequate to fully characterize the dynamics of DC assembly, modelling these reactions enables the prediction of system's behavior under mutations. On the other hand, the most important lack of the model is that it does not distinguish between the nucleus and cytoplasm. [START_REF] Mazemondet | Elucidating the Sources of β-Catenin Dynamics in Human Neural Progenitor Cells[END_REF] developed a model based on the [START_REF] Lee | The roles of apc and axin derived from experimental and theoretical analysis of the wnt pathway[END_REF] model that distinguishes the two compartments and also includes the β-catenin/Axin feedback loop. However, the specific model simplifies the DC assembly procedure. Here, we combine the two models.

The new model consists of all the chemical reactions in the [START_REF] Lee | The roles of apc and axin derived from experimental and theoretical analysis of the wnt pathway[END_REF] model except for R 13 (Fig. 1) which is replaced by the reaction

R 13 : β-cat n + T CF α19 ---- α20 β-cat n /T CF
where β-cat n corresponds to the β-cateninin molecules localized in the nucleus. The model is completed with one reaction for the β-cateninin transportation

R 14 : β-cat α23 ---- α24 β-cat n
and one for the regulation of the expression of Axin by the β-cat n /T CF complex

R 15 : β-cat n /T CF α25 β-cat n /T CF + Axin.
dx 1 dt = -α 1 x 1 + α 2 x 2 (16) dx 2 dt = α 1 x 1 -α 2 x 2 (17) dx 3 dt = α 3 x 4 -α 4 x 3 -α 10 x 11 x 3 + α 11 x 8 + α 13 x 9 (18) dx 4 dt = -α 5 x 2 x 4 -α 3 x 4 + α 4 x 3 + α 6 x 5 x 6 -α 7 x 4 (19) dx 5 dt = α 5 x 2 x 4 -α 6 x 5 x 6 + α 7 x 4 (
dx 15 dt = α 21 x 11 x 7 -α 22 x 15 (29)

Results

An implementation of the ABC-SMC method (Algorithm 2) is available in the Python package ABC-SysBio [START_REF] Liepe | ABC-SysBio-approximate Bayesian computation in Python with GPU support[END_REF]. Several kernels are already available in the ABC-SysBio. We used the component-wise Gaussian (cwG), the multivariate Gaussian (mvG) and the K-neighbors multivariate Gaussian (KmvG) for comparison. Assume a set of parameter vectors Θ Θ Θ = (θ 1 θ 1 θ 1 , . . . , θ N θ N θ N ) with weights w w w = (w 1 , . . . , w N ), where each parameter vector has d components θ i θ i θ i = (θ i 1 , . . . , θ i d ). In the cwG kernel each component θ j of the parameter vector is perturbed independently according to a Gaussian distribution with mean θ j and variance:

σ 2 j = N i=1 w i (θ i j -µ j ) 2 V -(V/V ) (33)
where V = N i=1 w i and V = N i=1 w i 2 are the sum of the weights and the sum of the squares of the weights respectively and µ j = N i=1 w (i) θ i j is the weighted mean for the specific component. In the mvG kernel, the correlations between the components of the parameter vector are taken into account. A particle θ j θ j θ j is perturbed using a multivariate normal distribution with mean θ j θ j θ j and covariance matrix:

Σ = V V 2 -V = N i=1 w i (θ i θ i θ i -µ µ µ) T (θ i θ i θ i -µ µ µ) (34)
where µ µ µ = N i=1 w i θ i θ i θ i is the weighted mean. The KmvG kernel, starts by selecting for each particle θ j θ j θ j the K-nearest neighbors. Then a multivariate normal distribution with mean θ j θ j θ j and covariance matrix computed by (34) but taking into account the K-nearest neighbors only, is used to draw a new particle. In ABC-SysBio, K is set by default to a quarter of the number of particles, and this value was used in the following studies. Further to the above described packages, featured by the ABC-SysBio framework, we have developed a Python package for the DPGMM and the Gibbs sampler described in Section 2. Regarding the DPGMM inference several iterations may be needed for the sampler to converge. In order to increase the overall computational efficiency, a reduced number of iterations (50) was pre-selected. At each iteration, the Bayesian information criterion (BIC) was used to evaluate how well the data fit to the proposed model. The cluster assignment with the lowest BIC was used to define the transition kernels in the ABC-SMC. For the ABC-SMC, we used 1000 particles and an Euclidean distance measure to compare simulated and experimental data. While different tolerance sequences may be optimal for each kernel, we used the same tolerance sequence for all of them as we wanted to compare their efficiency under the same conditions.

Toy example

In order to highlight the advantages of the DPGMM kernel, we initially apply ABC-SMC to a toy example with known posterior. We assume that the data are produced by the likelihood function:

p(D|θ θ θ) = N [θ 2 1 θ 2 2 ],
0.5 0 0 0.5 (35) and that D o = [1 1] is observed. The posterior density is then

p(θ 1 , θ 2 |D s ) = φ [1 1] θ 2 1 θ 2 2 , 0.5 0 0 0.5 (36)
where φ(•|µ µ µ, Σ) is the two dimensional normal density with mean µ and covariance matrix Σ. As it is shown in Figure 3A, there are 4 modes in the posterior distribution around the points

[1 1], [1 -1], [-1 1] and [-1 -1].
Using ABC-SMC and a uniform prior distribution on the space [-10 10]×[-10 10] we tried to estimate the posterior distribution. As it is shown in Figures 3C-F, the posteriors predicted by all 4 kernels are close to the true posterior. However, there is a significant difference in the efficiency of each kernel. When the cwG and mvG kernels are used many more particles have to be generated in order to accept 1000 particles from the next intermediate distribution in comparison to the locally concentrated KmvG and DPGMM kernels. The acceptance rates of the latter are in some cases over two times higher than the former ones. This is expected as the parameters of cwG and mvG kernels are derived from all the particles in the previous population and thus more and more of the generated particles fall into the area surrounded by the 4 modes. Moreover, the performance of the two kernels is almost identical, which is normal given that there is no correlation between the two parameters, due to the diagonal covariance matrix. From the other two kernels, DPGMM consistently provides better acceptance rates than KmvG. Especially in the last 5 populations, in which the acceptance rate is low for all kernels, 50 -70% more particles have to be generated with the KmvG than with the DPGMM kernel.

Parameter estimation for the Wnt pathway

The Lee-model (Section 3.1) consists of 22 parameters. While several of them have been measured through experiments on the Xenopus frog [START_REF] Lee | The roles of apc and axin derived from experimental and theoretical analysis of the wnt pathway[END_REF], we expect their values to be different for measurements performed on Human cells [START_REF] Tan | Analysis of wnt signaling β-catenin spatial dynamics in hek293t cells[END_REF]. Nevertheless, given the scarce experimental data available, we fixed all parameters to constant values that have been used in previous studies [START_REF] Maclean | Mathematical and statistical techniques for systems medicine: The wnt signaling pathway as a case study[END_REF] except for three: the binding rate between β-catenin and the destruction complex (α 11 ), the β-catenin degradation rate (α 16 ), and the binding rate of β-catenin and TCF (α 19 ). Uniform priors with ranges [0 100], [0 1], [0 1.5] were used for each parameter respectively. We also aimed at estimating the initial concentration of β-catenin. For the remaining species we used the initial concentrations reported in [START_REF] Tan | Wnt signalling pathway parameters for mammalian cells[END_REF]. For β-catenin, x 11 , we used a uniform prior with range [0.5 200]. A similar approach was used for the extended Lee model (Section 3.2). For the new parameters introduced due to the new reactions, the parameter values reported in [START_REF] Mazemondet | Elucidating the Sources of β-Catenin Dynamics in Human Neural Progenitor Cells[END_REF] were used. As with the Lee model all parameters were fixed except for: the binding rate between β-catenin and the destruction complex (α 11 ), the β-catenin degradation rate (α 16 ), and the β-catenin trasportation rate α 23 . For the first two the same uniforms as before were used, while for α 23 we used a uniform prior with range [0 10]. Given that two pools of β-catenin are considered in the extended model, we estimated the initial concentrations of both cytosolic and nuclear β-catenin, x 11 and x 16 . The simulated data were compared to the experimental data from [START_REF] Tan | Analysis of wnt signaling β-catenin spatial dynamics in hek293t cells[END_REF]. In [START_REF] Tan | Analysis of wnt signaling β-catenin spatial dynamics in hek293t cells[END_REF] the β-catenin concentration in the cytosol and the nucleus was measured at four time points (0, 60, 120, and 240 minutes) with and without Wnt. Here we consider the total β-catenin concentration after pathway activation. For the Lee model simulated data correspond to the concentration x 11 , while for the extended model we used the total β-catenin concentration, i.e. x 11 + x 16 . 35). (B) Acceptance rates (fraction of generated particles that are accepted) for the cwG purple), the mvG (cyan), the KmvG (green) and the DPGMM (red) kernels for each population. (C-F) Intermediate distributions for different tolerance values, starting with a high tolerance (red) to a low tolerance (purple) approximated using cwG (C), mvG (D), KmvG (E) and DPGMM (F) kernels.

Being mostly interested on the results than comparing the different kernels, we did not pre-define a tolerance sequence, but we used an adaptive approach to estimate the next tolerance based on the particles accepted at the current population. More precisely, after having accepted N particles for which the distance between simulated and experimental data is lower than the current tolerance i , we select the next tolerance i+1 to be equal to the higher distance of the αN best particles, where 0 < α < 1. This means that the αN particles from the current population will also be in the next population. A value of α close to 1 guarantees that the intermediate target distributions do not differ significantly. We set the α equal to 0.75 and the lowest tolerance T = 0.15. Based on the results from the toy example we used the DPGMM and the K-nearest neighbors kernels to ensure that the results will be the same for the two kernels.

The posterior distributions, as approximated using the DPGMM kernel, are shown in Figures 4 and5 for the model parameters and the concentrations respectively. Similar results were obtained using the K-nearest neighbors kernel. For the β-catenin-DC binding rate (α 11 ) the results are different for the two models. For the [START_REF] Lee | The roles of apc and axin derived from experimental and theoretical analysis of the wnt pathway[END_REF] model it spans the prior range, while for the compartmentalized model it takes values over the higher half of its prior range. The β-catenin degradation rate (a 16 ) values are concentrated towards the lower half for both models but the range for the Lee model is much wider. Finally, the values of the third parameter, the binding rate of β-catenin and TCF (α 19 ), span the prior range, while the β-catenin trasportation rate α 23 is concentrated around the upper limit of the prior range. It is clear from the results that some parameters affect more the results than others, for which their values do not signifficantly alter the results (i.e. the parameters whose values span the whole prior range). In summary for both models there exist sets of parameters for which model predictions match the experimental data well even if the remaining parameters have been fixed most probably to unrealistic values. This is not unreasonable if we consider the complexity of the models with respect to much more information than the parameter estimations. An interesting observation is that the initial β-catenin concentration for the Lee model has to be very low (the posterior is concentrated around 1, Fig. 5A) for the model predictions to be close to the experimental data. The posterior distributions for the compartmentalized model are significantly more spread. The concentration of the cytosolic β-catenin ranges from 0 to 50 nM, with a small peak around 20-25 nM, while the one for the β-catenin in the nucleus varies between 70 and 130 nM with a peak around 120 nM. These numbers are in accordance with the initial concentrations (25 nM in the cytosol and 147 nM in the nucleus) that have been exerimentally measured in HEK293T cells [START_REF] Tan | Wnt signalling pathway parameters for mammalian cells[END_REF].

Regarding the performance of the two kernels, the proposed DPGMM kernel provides either better or similar acceptance rates with the KmvG kernel (Fig 6A). However, the running times of the ABC-SMC algorithm (Fig. 6B) with the KmvG kernel are better than DPGMM. The reason is that the implemented Gibbs sampler is computationally expensive. However, it is important to note that the computational cost of simulating the data was extremely low for the model we consider here. One simulation of the ODE model lasted around 0.042 seconds. Stiffer or stochatic models usually require considerably more time. Projections of the running times of the ABC-SMC algorithm with acceptance rates equal to the ones observed using DPGMM and KmvG kernels and assuming simulations lasting 1 sec and 5 secs are

B A

Figure 6: Performance analysis using experimental data. (A) Acceptance rates (fraction of generated particles that are accepted) using the KmvG (red) and the DPGMM (black) kernels (B) Computational times for the two kernels (solid lines). Projection of the computational times assuming a simulation time of 1 sec (dashed lines) and 5 secs (dotted lines).

shown in Figure 6B (dashed and dotted lines respectively).

Conclusions

We presented a likelihood-free methodology to learn the parameters of mathematical models of complex biological systems. The proposed approach is based on Dirichlet process mixtures (DPMs) that have been previously used in density estimation. We used the DPMs in a ABC-SMC framework to efficiently explore the parameter space. It was shown that the proposed approach can bring potential benefits compared to commonly used alternatives improving the efficiency of SMC mathods in targeting complex, multimodal distributions. Apart from the improved performance, another advantage of DPMs is that they are user friendly as, due to their non-parametric nature, they do not require a predefined number of Gaussian components. The advantages of the method were highlighted using both a toy example and a real biological system, the Wnt signalling pathway.

While the proposed methodology was used with ABC-SMC its applicability is not limited to the specific framework. It can be easily extended to several bayesian inference algorithms, like MCMC, to design proposal distributions. Similarly, it could be used for other problems like model selection. In the ABC-SysBio package an ABC-SMC algorithm for model selection is available. Thus, the use of DPMs for model selection is straight forward. An interesting extension of the present work is the development of an ABC-SMC algorithm to estimate model parameters using single-cell experimental data. We envisage that the increased efficiency of DPMs in estimating multimodal distributions will be beneficial in studying cell-to-cell variability and infering multidimensional parameter distributions describing a population of cells that may be grouped according to their environmental conditions.

A drawback with DPM kernels is that they may become computationally expensive. The implemented Gibbs sampler requires a significant amount of time to infer the DPGMM. However, it is important to note that the computational cost of simulating biological systems may also be high. Thus, the increase in the acceptance rate may be beneficial when the computational cost to simulate complex data is not negligible. Apart from Gibbs sampling, several alternative inference approaches for DPMs are available [START_REF] Blei | Variational inference for Dirichlet process mixtures[END_REF][START_REF] Wang | Fast Bayesian inference in Dirichlet process mixture models[END_REF] and it would worth investigating their efficiency.

  the sum of the importance weights of the particles C j-1 c = {k = 1, . . . , N |z j-1 k = c} that belong to the component.

Figure 1 :

 1 Figure 1: Chemical reactions in the Lee et al. (2003) model.

Figure 2 :

 2 Figure 2: Differential equations in the Lee et al. (2003) model.

Figure 3 :

 3 Figure 3: Parameter estimation results for the toy example. (A) True posterior distribution for the likelihood function in (35). (B) Acceptance rates (fraction of generated particles that are accepted) for the cwG purple), the mvG (cyan), the KmvG (green) and the DPGMM (red) kernels for each population. (C-F) Intermediate distributions for different tolerance values, starting with a high tolerance (red) to a low tolerance (purple) approximated using cwG (C), mvG (D), KmvG (E) and DPGMM (F) kernels.

Figure 4 :Figure 5 :

 45 Figure 4: Marginal posterior distributions for the Lee et al. (2003) model parameters (A-C) and the compartmentalized model (D-F) approximated using ABC-SMC with DPGMM transition kernels.

Table 1 :

 1 [START_REF] Mazemondet | Elucidating the Sources of β-Catenin Dynamics in Human Neural Progenitor Cells[END_REF] the negative feedback loop formed through the regulation of Axin expression by nuclear β-catenin is modeled explicitly. In[START_REF] Jensen | A wnt oscillator model for somitogenesis[END_REF] the Wnt induced oscillatory Core species of Lee Wnt/β-catenin pathway model

			Species
		Acronym	var. ID	description
		Dsh a	x 1	active dishevelled
	Monomers	Dsh i GSK3 APC β-cat * β-cat	x 2 x 5 x 7 x 10 x 11	inactive dishevelled glycocen synthase kinase adenomatous polyposis coli phosporylated beta catenin beta catenin
		Axin	x 12	axin protein
		TCF	x 13	T-Cell factor
		APC * /Axin * /GSK3	x 3	complex
	Polymers	APC/Axin/GSK3 APC/Axin APC * /Axin * /GSK3/β-cat APC * /Axin * /GSK3/β-cat *	x 4 x 6 x 8 x 9	complex complex complex complex
		TCF/β-cat	x 14	complex
		APC/β-cat	x 15	complex

dynamics of a β-catenin/Axin negative feedback loop have been unraveled in ODEs models of somatogenesis. Most of the existing Wnt pathway models are based on the seminal work by

[START_REF] Lee | The roles of apc and axin derived from experimental and theoretical analysis of the wnt pathway[END_REF]

, who proposed an ODE, mono-compartment, core model of canonical Wnt pathway, validated through in-vitro measurement on the eggs of the Xenopus frog. Because of the amount of wet-lab experimental data, to date, the Lee model represents the most accurate/complete effort in Wnt pathway modelling.

  The differential equations for Axin and β-catenin become dx 11 dt = -α 10 x 11 x 3 + α 11 x 8 + α 15 -α 16 x 11 -α 21 x 11 x 7 +α 22 x 15 -α 23 x 1 1 + α 24 x 16

		(30)
	dx 12 dt	= -α

8 x 7 x 12 + α 9 x 6 + α 25 x 16 -α 18 x 12 (31) where x 16 denotes the concentration of β-catenin in the nucleus governed by the following differential equation: dx 16 dt = -α 19 * x 16 * x 13 + α 20 * x 14 + α 23 x 11 -α 24 x 16 (32)

  20) dx 6 dt = α 5 x 2 x 4 -α 6 x 5 x 6 + α 7 x 4 + α 8 x 7 x 12 -α 9 x 6 10 x 11 x 3 + α 11 x 8 + α 15 -α 16 x 11 -α 19 x 11 x 13 + α 20 x 14 -α 21 x 11 x 7 + α 22 x 15 (25) dx 12 dt = -α 8 x 7 x 12 + α 9 x 6 + α 17 -α 18 x 12 (26) dx 13 dt = -α 19 x 11 x 13 + α 20 x 14 -α 21 x 11 x 7 + α 22 x 15 (27) dx 14 dt = α 19 x 11 x 13 -α 20 x 14 (28)

			(21)
	dx 8 dt	= α 10 x 11 x 3 -α 11 x 8 -α 12 x 8	(22)
	dx 9 dt	= α 12 x 8 -α 13 x 9	(23)
	dx 10 dt	= α 13 x 9 -α 14 x 10	(24)
	dx 11 dt	= -α