Search tips
Search criteria 


Logo of ploscompComputational BiologyView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
PLoS Comput Biol. 2017 July; 13(7): e1005688.
Published online 2017 July 28. doi:  10.1371/journal.pcbi.1005688
PMCID: PMC5550005

A hierarchical Bayesian model for understanding the spatiotemporal dynamics of the intestinal epithelium

Oliver J. Maclaren, Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing,1,* Aimée Parker, Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft,2 Carmen Pin, Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing,2 Simon R. Carding, Funding acquisition, Project administration, Resources, Supervision,2,3 Alastair J. M. Watson, Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Writing – review & editing,2,3 Alexander G. Fletcher, Conceptualization, Methodology, Supervision, Writing – review & editing,4,5 Helen M. Byrne, Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing,6 and Philip K. Maini, Conceptualization, Funding acquisition, Project administration, Supervision, Writing – review & editing6
Qing Nie, Editor


Our work addresses two key challenges, one biological and one methodological. First, we aim to understand how proliferation and cell migration rates in the intestinal epithelium are related under healthy, damaged (Ara-C treated) and recovering conditions, and how these relations can be used to identify mechanisms of repair and regeneration. We analyse new data, presented in more detail in a companion paper, in which BrdU/IdU cell-labelling experiments were performed under these respective conditions. Second, in considering how to more rigorously process these data and interpret them using mathematical models, we use a probabilistic, hierarchical approach. This provides a best-practice approach for systematically modelling and understanding the uncertainties that can otherwise undermine the generation of reliable conclusions—uncertainties in experimental measurement and treatment, difficult-to-compare mathematical models of underlying mechanisms, and unknown or unobserved parameters. Both spatially discrete and continuous mechanistic models are considered and related via hierarchical conditional probability assumptions. We perform model checks on both in-sample and out-of-sample datasets and use them to show how to test possible model improvements and assess the robustness of our conclusions. We conclude, for the present set of experiments, that a primarily proliferation-driven model suffices to predict labelled cell dynamics over most time-scales.

Author summary

The intestinal epithelium is an important model system for studying the dynamics and regulation of multicellular populations. It is characterised by rapid rates of self-renewal and repair; dysregulation of these processes is thought to explain, in part, why many tumours form in the intestinal and similar epithelial tissues. These features have led to a large amount of work on estimating cell kinetic parameters in the intestine. There remain, however, large gaps between the raw data collected, the interpretation of these experimental data, and mechanistic models that describe the underlying processes. Hierarchical statistical modelling provides a natural method with which to bridge these gaps, but has, to date, been underutilised in the study of intestinal tissue self-renewal. As we illustrate, this approach makes essential use of the distinction between ‘measurement’, ‘process’ and ‘parameter’ models, giving an explicit framework for combining experimental data and mechanistic modelling in the presence of multiple sources of uncertainty. We apply this approach to analyse experiments on healthy, damaged and recovering intestinal tissue, finding that observed data can be explained by a model in which cell movement is driven primarily by proliferation.



The intestinal epithelium provides crucial barrier, transport and homeostatic functions. These requirements lead it to undergo constant repair and regeneration, and dysfunctions can result in pathologies such as tumorigenesis [17]. Much work has been carried out on estimating the rate parameters in the intestine and other epithelia [1, 810]. However, attempts to interpret these experimental data using mechanistic modelling remain inconclusive [1114], due to the lack of consistent and reproducible approaches for comparing models representing conjectured biological mechanisms, both to each other and to experimental data.

This challenge goes in both directions: using experimental data (taken to be ‘true’) to parameterise and test mathematical or computational formalisations of mechanistic theories, and using these models (taken to be ‘true’) to predict, interpret and question experimental results. Both experimental measurements and mathematical models are subject to uncertainty, and we hence need systematic ways of quantifying these uncertainties and attributing them to the appropriate sources. Furthermore, establishing a common approach for analysing experimental results, formulating mechanistic models and generating new predictions has many potential advantages for enabling interdisciplinary teams to communicate in a common language so that they may efficiently discover and follow promising directions as and when they arise.


We address the above issues by developing a hierarchical Bayesian model for combining measurements, models and inference procedures, and applying it to a set of experiments targeting mechanisms of repair and regeneration in the intestinal epithelium.

While progress is now being made in tackling this challenge in other areas of biology (see, for example, [1517]), to our knowledge the problem of intestinal epithelial dynamics has not yet been investigated using such an approach.

The experiments under investigation were performed by ourselves and are presented in more detail in [18]. The aim of these experiments was to determine how proliferation rates, tissue growth and cellular migration rates are related under healthy, damaged (Ara-C treated) and recovering conditions, and how these relations can be used to identify mechanisms of repair and regeneration.

A notable feature of the Bayesian approach to probabilistic modelling is that all sources of uncertainty are represented via probability distributions [1921]. Adopting this perspective, we consider both observations and parameters to be random variables. Within a modelling context, uncertainty may be associated with (at least): parameters within a mechanistic model of a biological or physical process, the mechanistic model of the process itself and the measurements of the underlying process. This leads, initially, to the postulation of a full joint probability distribution for observable, unobservable/unobserved variables, parameters and data.

Another key feature of the Bayesian perspective is that it provides a natural way of decomposing such full joint models in a hierarchical manner, e.g. by separating out processes occurring on different scales and at different analysis stages. A given set of hierarchical assumptions corresponds to assuming a factorisation of the full joint distribution mentioned above, and gives a more interpretable and tractable starting point.

Our factorisation follows that described in [2124]. This comprises: a ‘measurement model’, which defines the observable (sample) features to be considered reproducible and the precision with which they are reproducible (the measurement scale); an underlying ‘process’ model, which captures the key mechanistic hypotheses of spatiotemporal evolution, and a prior parameter model which defines a broad class of a priori acceptable parameter values.

In order to illustrate some of the modelling benefits of the hierarchical approach, we show how both discrete and continuous process models can be derived and related using the hierarchical perspective. We discuss the relationship between the conditional/hierarchical modelling and causal modelling literatures (see [2527] for reviews) and illustrate the distinct roles of (Bayesian) predictive distributions vs. parameter distributions for model checking and the assessment of evidence, respectively [20, 2831].


Our hierarchical Bayesian model incorporates measurement, process and parameter models and facilitates separate checking of these components. This checking indicates, in the present application to the spatiotemporal dynamics of the intestinal epithelium, that much of the observed measurement variability can be adequately captured by a simple measurement model. Similarly we find that a relatively simple process model can account for the main spatiotemporal dynamics of interest; however, model checking also identifies a minor misfit in the process model appearing over long time-scales. This motivates possible model improvements: we consider the inclusion of additional finite-cell-size effects in the process model, derived from a discrete process model and a subsequent continuum approximation formulated in terms of conditional probability. This only gives a slightly better qualitative fit to experimental data, however. We instead find that the dominant sources of the long-time misfits are probably due to some other factors—most likely relatively slow, time-varying proliferation rates (e.g. due to circadian rhythms). In summary, a primarily proliferation-driven model appears adequate for predictions over moderate time-scales.

Materials and methods

Experimental treatments and data processing

Homeostasis mouse model

To estimate intestinal epithelial cell proliferation and migration rates under normal, homeostatic conditions in healthy mice, we used standard methods of proliferative cell labelling and tracing [1, 810, 3234] (see also [18] for full details). Actively proliferating cells in the intestinal crypts were labelled by single injection of the thymine analogue 5-bromo-2-deoxyuridine (BrdU) and labelled cells detected by immunostaining of intestinal sections collected from different individuals over time. Migration of labelled cells traced from the base of crypts to villus tips was monitored over the course of 96 hours (5760 min). At least 30 such strips were analysed per mouse. The figures presented in [18] reveal that these strips were independent and obtained from one-cell thick sections. All strips in which the base of the crypt and the tip of the villus were clearly observed were considered. All sides of the crypt that were visible and connected to entire villi were analysed. A typical image from those analysed in [18] is shown in Fig D in S1 Supplementary information.

Blocked-proliferation mouse model

To assess the effects of proliferation inhibition on crypt/villus migration, migrating and proliferating epithelial cells were monitored by double labelling with two thymine analogues (BrdU and IdU), administered sequentially a number of hours apart and subsequently distinguished by specific immunostaining in longitudinal sections of small intestine. Following initial IdU labelling of proliferating cells at t = −17h (−1020 min, relative to Ara-C treatment), mice were then treated with cytosine arabinoside (Ara-C) at 250 mg/kg body weight, a dose reported to inhibit proliferation without causing major crypt cell atrophy (see [18] and references therein for details). Tissues were collected over 24 hours, with BrdU administered one hour prior to collection to check for residual proliferation. Successful inhibition of proliferation following treatment with Ara-C was confirmed by an absence of BrdU (S-phase) and phospho-Histone H3 (pH3) staining (M-phase) in longitudinal sections of small intestine (again, see [18] for full details).

Recovering-proliferation mouse model

The above Ara-C treatment effect was observed to last for at least 10h (600 min). Cell proliferation returned to near normal levels in samples obtained 18h (1080 min) post-Ara-C treatment. We hence considered samples collected at least 10h post-Ara-C treatment as corresponding to ‘recovering-proliferation’ conditions.

Data processing: Reference grid and key observable features

To relate experimental measurements to the mechanistic models discussed below we specified a reference grid and defined the key features of the data relative to this grid. These key features established an ideal ‘underlying population’ from which samples were considered to be drawn. This also allowed us to construct our hierarchical model in a (data-to-parameter) manner, starting from a measurement model.

With reference to Fig 1, we viewed the data as a collection of one-dimensional ‘strips’ of cells. The strips extended from the base of the crypt to the tip of the villus, along the so-called ‘crypt-villus’ axis. This corresponds to how strips were collected experimentally, but does not account for possible biases due to ‘angled’ sampling [35, 36]. Each measurement was given a spatial cell location index i and a time label t. The location index was measured in numbers of cells along the crypt-villus axis, starting from the crypt base, and hence defined a discrete one-dimensional grid.

Fig 1
The (a) intestinal epithelium, (b) individual measurements as strips of cells and (c) collection of strips, where ‘C’ and ‘V’ indicated ‘crypt’ and ‘villus’ respectively.

When notationally convenient, the labels i and t were combined in a two-dimensional grid of space-time points via the single index parameter s := (i, t). A ‘typical’ reference crypt-villus unit was characterised by the two vectors (L, n), where L is the vector of underlying labelled fractions (i.e. occupancy probabilities) at each grid point and n is the vector of number of samples at each grid point. This defines a useful reduction of the system from two spatial dimensions to one.

We assumed that each strip was independent of the others as, in general, strips were taken from different crypt-villus units and/or animals after ‘identical preparation’. Thus we did not ever directly possess, for example, measurements of a particular crypt with dimensions given in terms of a certain number of strips. We note, however, that the dynamics of strips in a given crypt may be affected by those in the same crypt. We did not consider this additional complexity in the present work, and so this complication should be kept in mind when interpreting the results.

Mathematical model

Our hierarchical probability model was constructed on the basis of conditional probability assumptions. These allowed us to factor out a measurement model, a mechanistic model and a parameter model.

Overall hierarchical structure

Our model structure consisted of a full joint distribution, conditioned on a given experimental treatment E and known sample size vector n, decomposed according to

p(yLk|nE) = p(y|Ln)p(L|k)p(k|E), 

where k are the cellular proliferation rates (these are discussed below). This hierarchical factorisation corresponds to the assumption of conditional independence between the various levels, i.e. p(y|L, k, n, E) = p(y|L, n), p(L|k, n, E) = p(L|k) and p(k|n, E) = p(k|E). The first term, p(y|L, n) is the measurement model; the second term p(L|k) is the underlying process model, and the last term p(k|E) is a prior parameter model. These are discussed in turn below.

Notably, a ‘causal’ (structural invariance) assumption [2527, 3742] is made by assuming that the experimental treatment condition affects the process parameters k but not the structure of the measurement or process models. In particular, while the experimental treatment ultimately affects y, in our model it does so via its effect on k and k’s subsequent effect on L. This leads to the conditional independence structure mentioned. In terms of so-called directed acyclic graphs (DAGs), used in the causal modelling literature referenced above, we assumed

E → k → L → y.

Another way of stating this is that knowledge of L (and n) is sufficient to determine y regardless of how L was brought about, but that to know L we (ultimately) have to know which experiment was carried out. Note that in general we suppressed, in our notation, the explicit conditioning on sample size n, since it was taken to be fixed and known, as well as the conditioning on E (keeping in mind that it was assumed to only affect k).

The assumptions underlying the above factorisation could be checked to some extent. This relied on a distinction between working ‘within’ the model—e.g. parameter estimation assuming the model and factorisation is valid—and working ‘outside’ the model, e.g. checking the validity of the model structural assumptions themselves [20, 30, 31]. This distinction is made in the Results section.

Implicit in the model derivations, discussed below, we used a deterministic expression of conservation of probability for the process model, as is typical for such equations [43]. It sufficed for the presentation here to simply replace all functional dependencies on the process variable above with a dependence on the process parameters [21].

We provide further discussion of our approach to interpretation of statistical evidence in the Supplementary Information.

Bayesian framework for predictions and incorporating information from observations

The overall model of the previous section defined our initial ‘generative’ probabilistic model, prior to explicitly incorporating information from our experimental data. This enabled samples to be drawn from both prior predictive and prior parameter models, in the usual way (see e.g. [20, 44] and the Computational methods section below). In particular, the prior predictive distribution was used in its usual form

y ∼ p(y) = ∫p(y|f(k))p(k)dk

which incorporates the aforementioned deterministic link between a given sample of process parameters and the output process variable, L = f(k). Note that here ~ denotes ‘distributed as’, or more relevantly, ‘samples drawn according to’.

To incorporate new data y0 we updated the parameters of the model, hence passing to a ‘posterior predictive’ model [20]

y|y0 ∼ p(y|y0) = ∫p(y|f(k))p(k|y0)dk

where we used the conditional probability closure assumption p(y|f(k), y0) = p(y|f(k)). This closure assumption can be interpreted as maintaining our same mechanistic model despite new observations.

Measurement model

The measurement model p(y|L, n) component was taken to be a binomial distribution of the form


This related our ‘raw’ observable y, the vector of counts of labelled cells at each grid point, to ‘ideal characteristics’ of comparison (L, n). We took the measurement component to be independent of the experimental treatment E, i.e. treatment was assumed to affect the underlying process parameters only (see ‘Overall hierarchical structure’ section). The measurement component p(y|L, n) then defined a likelihood function for this measurement model,


When interpreting model misfit based on residuals, we applied the usual normal approximation to the binomial distribution. In that case, denoting the set of all measured labelled fractions through the (useful, but slightly non-standard) notation y/n := (y1/n1, …, yS/nS), we have


where the standard deviations are given by σs=Ls(1-Ls)ns.

Process models: Spatially discrete model

We developed our process model at different levels of resolution. First, we considered a spatially discrete probabilistic model at the level of our measurement grid defined above. Second, we considered two spatially continuous approximations: one model excluding explicit cell-scale effects and one including them. Here we first consider the discrete model.

Our basic ‘process’ model described the evolution of the occupancy probabilities (population labelled fractions) at the scale of the measurement grid. This was derived as follows. With reference to Fig 1, we considered a collection of one-dimensional ‘strips’ of cells. We used li [set membership] {0, 1} as an indicator variable denoting the occupancy status of site i of a given strip. The full state of this strip was given by the vector l = (l0, l1, …lS − 1). We then sought a description of the probabilistic dynamics in terms of a discrete-time Markov chain for the probability distribution of the full state p(l, t) following standard arguments [43, 45].

We began from an explicit joint distribution for the full state and then reduced it to a description in terms of the set of ‘single-site’ probability distributions p(li, t) for each site i. This derivation was aided by adopting an explicit notation: the probabilities of occupancy and vacancy at site i at time t were denoted by p(li(t) = 1) and p(li(t) = 0) respectively. Since p(li(t) = 1) + p(li(t) = 0) = 1 we only needed to consider the probability of occupancy to fully characterise the distribution p(li(t)).

The equation of evolution for this probability was derived by considering conservation of probability in terms of probability fluxes in and out, giving, to first order in Δt,


The first term on the right represents a net ‘influx of occupancy probability’ due to a single division event at site j < i, each division event having a probability kj Δt. This flux means the value of the state variable li(t) = 0 could be replaced, at the next time step, by the value of li − 1(t) = 1. Similarly the second term represents a net ‘outflux of occupancy probability’ due to a division event at site j < i. Partitioning on the events li − 1(t) = 0 and li(t), we obtain


Process models: Underlying continuous model and zeroth-order approximation

To aid model interpretation and model cross comparisons we derived a continuous approximation to the occupancy probability, L(x, t). This gave a further idealisation of the ‘underlying population’ from which we envisaged the strips were sampled. This smoothness assumption, while not strictly necessary, meant some model properties could be interpreted in terms of local derivatives; it also reduced arbitrary dependence on discrete grid features, aiding future comparisons with off-lattice and/or continuum models [44].

As detailed in the Supplementary Information, L(x, t) satisfies the equation




where k(x) denotes the (net) proliferation rate.

Process models: Underlying continuous model and higher-order spatial effects

Our ‘zeroth-order’ continuous approximation above was obtained by neglecting all higher-order terms in Δx. We anticipate that a more accurate continuum approximation may be obtained by retaining higher-order spatial derivatives and hence finite-cell-size effects. This gives rise to a Fokker-Planck equation containing a diffusion term [43]:


where D(x) = (1/2)Δxv(x). Similar equations have been derived before, based on continuous approximations to discrete master equations [4649]. Retaining the second spatial derivative hence amounts to accounting for spatial effects due to finite cell sizes. We evaluated our original ‘zeroth-order’ (advection) model against our data, and also examined the extent to which higher-order spatial terms such as those considered above could account for any misfits.

Definition of priors

Since we adopted a Bayesian perspective in this work we required a parameter prior model that could express additional modelling assumptions [20]. Essentially, this is an empirical Bayesian approach as we used an empirical correlation matrix [50].

Candidate proliferation profiles, varying with cell locations, were represented as realisations from a prior given in terms of a discretised random field (a random vector) k of length m = 5, modelled as a multivariate Gaussian 𝒩(μ, C) with joint distribution


characterised by its mean vector μ and covariance matrix C. This parameter prior constrained the variability of the spatially varying parameter field a priori to help avoid unphysical solutions.

The covariance matrix was first decomposed into a standard deviation matrix given by the outer (tensor) product of the standard deviation vector for each variable, S = σσT, and correlation matrix R. These multiply element-wise to give Cij = Sij Rij (no summation). We then adopted the common, equivalent, representation C = DRD where D is a diagonal matrix with diagonal entries Dii = σi.

We took the correlation matrix R to have the squared-exponential (Gaussian) correlation function k(i,j)=exp((i-j)22lc2), where lc is a parameter controlling the characteristic length-scale of the correlations in terms of number of indices of k. This characteristic length scale gives the number of k indices over which the correlation function decays to 1/e. This allowed us to control the ‘smoothness’ of the realisations from the k prior, in the sense that as lc is increased the values ki and kj tend to be more similar.

The matrix R was generated by evaluating this correlation function at discrete locations along the crypt-villus axis. This discretisation was chosen to be coarser than the measurement grid and gave a variation somewhat similar to compartment-style regions of proliferation activity. This corresponded to assuming that the cell-type and associated proliferation rates varied on a coarser scale than individual cells, and was thus somewhat similar to a compartment-style assumption [51, 52], though the resulting proliferation rate function is defined for all values of the finer, individual-cell scale x. The parameter lc could also be interpreted as a ‘parameter correlation length’ for the proliferation rates, a measure of the number of parameters—or number of ‘compartments’—over which the correlations decay. We considered correlation lengths of 1–2 parameters.

We found it most informative to visualise realisations of the whole function from the resulting prior rather than simply give the individual parameters/matrices separately ([21] discusses this visualisation approach to priors in more detail). These are hence discussed and displayed in more detail in the Results section below.

Implementation of MCMC sampling and Bayesian updating

We used the (open source) Python package emcee ( to perform Markov Chain Monte Carlo (MCMC) [53] to obtain samples from the posterior distributions. Given samples from the resulting prior and posterior parameter distributions, respectively, prior and posterior predictive distributions were obtained by forward simulation of the process model described below. We note that each candidate proliferation rate vector k is connected to the measurements y via the latent vector L; since this step is deterministic, however, no additional sampling steps were required for the process model component.

Solution of differential equations

For the results in all sections other than the final results section in which we include higher-order spatial effects, we solved the differential equation model using the PyCLAW [54, 55] Python interface to the CLAWPACK [56] set of solvers for hyperbolic PDEs. We adapted a Riemann solver for the colour equation available from the Riemann solver repository ( For testing the inclusion of higher-order spatial effects (thus changing the class of our equations from hyperbolic to parabolic) we used the Python finite-volume solver FiPy [57].

Data and source code availability

Our code is available in the form of a Jupyter Notebook ( in the Supplementary Information. We ran these using the Anaconda distribution of Python ( which is a (free) distribution bundling a number of scientific Python tools. Any additional Python packages and instructions which may be required are listed at the beginning of our Jupyter Notebook.


Parameter inference under homeostatic (healthy) conditions

Fig 2 illustrates the process of updating from realisations of the prior distributions of the proliferation and velocity fields to realisations of their posterior (post-data) distributions. As discussed in the Materials and methods section above, these are generated by an underlying piecewise-constant Gaussian random field of proliferation rates, k. This has length m = 5, and defines an assignment of the cell indices into biologically-motivated regions of proliferation activity.

Fig 2
Simulated realisations from the prior (left) and posterior (right) distributions for proliferation profiles (top) and velocities (bottom).

The left-hand side of the figure shows simulations from the prior distribution for proliferation field (top) and realisations from the induced distribution for the velocity field (bottom), respectively. The right-hand side shows the corresponding simulations after the prior parameter distribution has been updated to a posterior parameter distribution. The prior-to-posterior parameter estimation was carried out using the MCMC sampling approach described above with t = 120 min (2 h) as an initial condition and t = 360 min (6 h) and 600 min (10 h) as given data. The initial condition for the underlying labelled fraction (occupancy probability) was determined by fitting a smoothing spline to the data. The prior distribution for the proliferation field shown in Fig 2 incorporated a weak mean trend in net proliferation rates, rising from the crypt base to the mid-crypt before falling exponentially to zero over the last few parameter regions post-crypt end, and a parameter correlation length of 1. These assumptions can be relaxed/varied with little effect, though typically a non-zero parameter correlation length and a shut-off in proliferation after the crypt end produce more stable (well-identified) estimates. Additional visualisations of the parameter inferences are provided in Fig A-C in S1 Supplementary information.

Parameter inference for blocked proliferation conditions

Fig 3 is the same as Fig 2 described in the previous section, but this time under treatment by Ara-C. Results from the baseline case are shown in grey, while those from Ara-C treatment are shown in blue. Here 1140 min (19 h post IdU labelling, 2 h post Ara-C treatment) was used as the initial condition and 1500 min (25 h post IdU labelling, 8 h post Ara-C treatment) used for fitting. The intermediate time 1260 min (21 h post IdU labelling, 4 h post Ara-C treatment) and later time 1620 min (27 h post IdU labelling, 10 h post Ara-C treatment) were used as out-of-sample comparisons (see later).

Fig 3
Simulated realisations from the prior (left) and posterior (right) distributions for proliferation profiles (top) and velocities (bottom) under Ara-C treatment (blue) as compared to no treatment (grey).

As can be seen, there is clear inhibition of proliferation and an even clearer effect on the cell migration (growth) velocity. The greater variability in the underlying parameter results compared to the baseline case may indicate, for example, greater parameter underdetermination and/or inconsistency of the model. This is not surprising as we expect all proliferation parameters to be reduced to similar (low) values and hence to become less distinguishable.

To add additional stability to the results we can attempt to reduce underdetermination in the parameters by increasing the parameter correlation length and inducing an effectively more ‘lumped’ representation of the parameter field (since values tend to stick together more). Doing this removed the more extreme negative net proliferation in the posterior profile, however it still allowed for small amounts of negative net proliferation/velocity (the available Jupyter notebook can be used to explore various prior assumptions).

Again, the need to introduce more stability is likely due to some combination of the limitations of resolution, a consequence of trying to fit the data too closely, or an indication of model inadequacies. In particular, under inhibited-proliferation conditions the effective number of parameters would be expected to be reduced. When fitting the full model, with largely independent parameters for each region, it is to be expected that some additional regularisation would be required for greater stability.

Parameter inference for recovering proliferation conditions

Ara-C is metabolised between 10–12 h post-treatment. The two times considered here, 1620 min and 2520 min, correspond to 10 h and 25 h post Ara-C treatment, respectively, i.e to the end of the effect and after the resumption of proliferation. Hence, to check for the recovery of proliferation, we fitted the model using 1620 min as the initial condition and 2520 min as the final time.

Fig 4 is the same as Figs Figs22 and and33 described in the previous sections, but this time after/during recovering from treatment by Ara-C. The previous results from the baseline case are shown in grey, while the new results following recovery from Ara-C treatment are shown in blue. Here 1620 min (27 h post IdU labelling, 10 h post Ara-C treatment) was used as the initial condition and 2520 min (42 h post IdU labelling, 25 h post Ara-C treatment) used for fitting. We did not make additional out-of-sample comparisons in this case, though in-sample posterior predictive checks were still carried out (see later).

Fig 4
Simulated realisations from the prior (left) and posterior (right) distributions for proliferation profiles (top) and velocities (bottom) after recovery from Ara-C treatment (blue) as compared to no treatment (grey).

Here, as expected, the proliferation and velocity profiles indicate that proliferation has resumed. The rates of proliferation appear to be lower than under fully healthy conditions, however, perhaps due to incomplete recovery (the initial condition being right at the beginning of the recovery period). The timing of the recovery of proliferation and the well-identified proliferation and velocity profiles inferred give no indication that any other mechanism is required to account for these data, however.

Predictive checks under homeostatic (healthy) conditions

Fig 5 illustrates simulations from the predictive distributions corresponding to the prior and posterior parameter distributions of Fig 2. This enables a first self-consistency check—i.e. can the model re-simulate data similar to that to which it was fitted [20, 58]? If this is the case then we can (provisionally) trust the parameter estimates in the previous figure; if not, then the parameter estimates would be unreliable, no matter how well-determined they seem. In our case the model appears to adequately replicate the data used for fitting.

Fig 5
Simulated realisations from prior (top) and posterior (bottom) predictive distributions (grey) for label data at fitted times (120 min, 360 min and 600 min i.e. 2 h, 6 h and 10 h).

Figs Figs66 and and77 illustrate two additional ways of visualising replicated datasets. The former visualises the label profile along the crypt-villus axis at the future unfitted/out-of-sample time 1080 min (18 h), while the latter visualises both fitted (120 min/2 h, 360 min/6 h and 600 min/10 h) and unfitted/out-of-sample (1080 min/18 h) predictions plotted in the characteristic plane (t, x) in which the slopes along lines of constant colour should be inversely proportional to the migration velocities at that point, due to the (hyperbolic) nature of our ‘colour’ equation (see e.g. [59]). We have interpolated between the dotted grid lines. These figures, in combination with Fig 5, indicate that the model is capable of reliably reproducing the data to which it was fitted, as well as predicting key features of unfitted datasets such as the rate of movement of the front. On the other hand, there is clearly a greater misfit with the predicted rather than fitted data. In order to locate the possible source of misfit we considered various model residuals and error terms—see ‘Locating model misfit’ below.

Fig 6
Simulated realisations from prior (left) and posterior (right) predictive distributions (grey) for label data at the unfitted (out-of-sample) time 1080 min (18 h).
Fig 7
Actual (smoothed) data (left, black box) and one replication based on the model (right; plotting the latent/measurement-error-free process) as visualised in the characteristic plane.

Predictive checks under blocked proliferation conditions

Here data from 1140 min (19 h; post IdU labelling) were used as the initial conditions and 1500 min (25 h) used for fitting. 1260 min (21 h) and 1620 min (27 h) were used as out-of-sample (non-fitted) comparisons. Fig 8 is analogous to Fig 5 in the healthy case. In general all of the features up to 1620 min (27 h) in Fig 8, and for both fitted and predicted times, are reasonably well captured. The fit at 1620 min is generally good, but perhaps worse than the other cases. This could be due to errors in longer-time predictions and to the beginning of proliferation recovery: we explore these alternatives in what follows.

Fig 8
Simulated realisations from posterior predictive distributions (grey) for label data at 1140 min (initial condition), 1500 min (fitted time) and at two out-of-sample/unfitted times (1260 and 1620 min).

Predictive checks under recovering proliferation conditions

As discussed above, Ara-C is metabolised between 10–12 h post-treatment. The two times considered here, 1620 min and 2520 min, correspond to 10 h and 25 h post Ara-C treatment, respectively, i.e to the end of the effect and after the resumption of proliferation.

As can be seen in Fig 9, and as expected, the label has resumed movement in concert with the resumption in proliferation. The model appears to fit reasonably well.

Fig 9
Simulated realisations from posterior predictive distributions (grey) for label data at 1620 min (initial condition) and 2520 min (fitted).

Locating model misfit

While the zeroth-order model behaves essentially as desired under experimental perturbation, and is likely capturing the essential features of interest, we observed some minor model misfit. We used posterior predictive checks to unpick the contributions of the various model parts and determine the source(s) of misfit. This in turn motivated potential model improvements. These checks were carried out under baseline (healthy) conditions as we were more confident of the experimental results under this scenario, but they can equally be carried out for the other datasets. Note, however, that time-varying effects are not expected to be as relevant under conditions of inhibited proliferation.

Fig 10 shows the following checks: measurement error as determined by subtracting a smoothed spline from the observed data (dark line) and comparing these to the results obtained by subtracting the process model from the simulated data (panels 1–4, moving left-to-right and top-to-bottom, showing fitted—120 min/2 h, 360 min/6 h and 600 min/10 h—and unfitted/out-of-sample—1080 min/18 h—times). This presentation follows the noise-checking approach in [60], as well as the general recommendations given in [20, 58]. Reliable interpretation of these as ‘true’ measurement residuals depends on the validity of the normal approximation Eq 7 since these expressions are not directly interpretable in terms of the discrete binomial model (see e.g. [20, 58]). These are also visualised in terms of the corresponding cumulative distributions in the middle panel (panel 5, following as above). Panels 6–9 show the differences between the underlying process model and the smoothed spline fitted to the data. As can be seen, the measurement model appears approximately valid at all times, while the process model appears to have non-zero error for the 1080 min sample. We consider this in more detail next.

Fig 10
Model and data residual components.

Possible model improvement and robustness—Higher-order spatial effects

As discussed in the process model section above, the presence of cellular structure in the epithelial tissue means that higher-order spatial effects could be present. One way of deciding whether these are important is to consider the extent to which these may account for the minor misfit identified above, as opposed to other factors such as time-varying proliferation rates. To do this we considered both uniform percentage reductions of the original parameter estimates (approximating time-varying rates) and the inclusion of higher-order spatial terms.

Fig 11 gives an idea of the qualitative differences induced by including the higher-order spatial terms and those that could be induced by time-varying proliferation rates. This figure is based on the (healthy) 1080 min (18 h) data in which we found some indication of a process model error.

Fig 11
Comparison of the modified process model which includes higher-order spatial terms (blue) to the original model (grey, dashed), both at lowered proliferation rates (decreased 20%), which is required for a better fit to the data.

We see that while the higher-order model appears to give a slightly better qualitative fit to the data, both the higher-order and lower-order models require similar reductions of the parameter values to quantitatively improve the fit to our out-of-sample data. The reduced parameter values shown in Fig 11 correspond to a reduction of 20%, which was chosen visually as a reduction accounting for the bulk of the misfit.

Thus the key (yet relatively small) difference between the model and out-of-sample data is likely due to an effect other than finite-cell sizes; in this case it is likely due to time-variation in parameter values due to circadian rhythms (we have assumed steady-state parameter values). Other potential factors include label dilution or an unmodelled mixing phenomenon in the full two-dimensional case. We note however that these effects are small and appear to be important primarily for predicting much further ahead in time than the fitted data and the steady-state parameter assumption is likely valid for reasonable time intervals. This means that the more easily interpretable original model may be sufficient for many purposes.


Understanding the intricate dynamics of the intestinal epithelium requires an interdisciplinary approach that integrates experimental measurements, mathematical and computational modelling, and statistical quantification of uncertainties. While a diverse range of mathematical models have been proposed for epithelial cell and tissue dynamics (reviewed in [51, 52, 6163]), from compartment models to individual-based models to continuum models, we lack consistent and reproducible approaches for comparing models representing conjectured biological mechanisms both to each other and to experimental data (for an overview, see our review [44]). These shortcomings may explain why questions such as the connection between proliferation and migration and its variation under experimental perturbations remain open [814].

The aim of the present work was to acknowledge and confront these difficulties explicitly, and to present some initial constructive steps towards establishing such a framework. To do this we carried out new experiments (described more fully in a companion paper [18]) aimed at determining how cell proliferation rates, tissue growth and cell migration rates are related in the intestinal epithelium under healthy, damaged (Ara-C treated) and recovering conditions. We performed BrdU/IdU cell-labelling experiments under these respective conditions. We then developed a probabilistic, hierarchical (conditional) model to process and interpret these data.

Our hierarchical model provides a best-practice approach for describing and understanding the uncertainties that could lead to unreliable mechanistic conclusions—uncertainties in experimental measurement and treatment, difficult-to-compare mathematical models of underlying mechanisms, and unknown or unobserved parameters. Our approach was influenced by recognising the benefits that the hierarchical Bayesian approach has demonstrated in applications across a number of different disciplines (e.g. in environmental and geophysical science as in [64, 65]; ecological modelling as in [66, 67]; and in Bayesian statistical modelling and inverse problems more generally as in [2024, 68]). We also note that a hierarchical approach can have significant benefits outside the Bayesian framework (see for example the ‘extended likelihood’ approach described in [6971]).

The hierarchical approach provides a framework, not only for combining disparate sources of uncertainty, but also for facilitating modelling derivations and relating discrete and continuous models. Though the resulting measurement, process and parameter models can be (or have been) derived by other means, as far as we are aware this particular perspective has not been systematically utilised in the manner considered here.

We also note the connection between the choice of a measurement model as required here (and/or process model error, and following e.g. [2124, 64, 72]), and the development of approximate sampling and parameter fitting procedures, which are particularly useful for analytically difficult models. A key concern of the latter is the appropriate choice of summary statistics for constructing a ‘synthetic likelihood’ [73] or similarly-modified posterior target for Approximate Bayesian Computation (ABC) [7476]. This choice determines (implicitly or explicitly) in which ways a given model or set of models can be considered an ‘adequate’ representation of the data, which features are considered to be reproducible and what the associated ‘noise’ structure should be ([77] presents an alternative approach to characterising data features and model adequacy). These issues are crucial in deciding how to model the complexity of epithelial cell and tissue dynamics.

An important next step, as described above, would be to consider other process model types and to evaluate and compare them under carefully modelled experimental conditions. Extensions incorporating other mechanical and/or cellular-level information (e.g. [11, 12]) into process models would provide a natural next step. Importantly, due to the separation between measurement and process model components, these more complex process models could be incorporated into our present framework simply by replacing our process model component with a new model, while retaining the same measurement model. Of course additional parameters would require additional prior assumptions, and if additional data features were of interest then these would need to be incorporated into a modified measurement model. The benefit of a hierarchical model is that it offers an explicit guide as to where such modifications should be incorporated.

In summary, the main results established using the above approach were

  • An adequate description of intestinal epithelial dynamics is achievable using a model based on proliferation-driven growth alone;
  • This model is consistent with healthy, proliferation-inhibited (Ara-C treated) and recovering conditions;
  • The measurement and process model errors can be reasonably distinguished and checked separately;
  • Checking indicates that much of the natural variability is attributable to the data collection process and this process can be modelled in a simple manner
  • Possible model errors can also be identified and proposed explanations incorporated and tested within our model, and, thus, the proper interpretation of experimental procedures is aided by using an explicit mathematical model and its predictive simulations
  • Including finite-cell-size effects gives a slightly better qualitative fit to experimental data, but the dominant sources of the long-time misfits are likely due to some other factor such as (relatively slowly) time-varying proliferation rates (e.g. due to circadian rhythms) or label dilution.

Supporting information

S1 Supplementary information

Fig A. Posterior for proliferation rates under baseline, healthy conditions. The upper diagonal represents the marginal distributions for each proliferation rate when averaging over all other profileration rates. The plots below the diagonal show bivariate marginal distributions illustrating pairwise associations after averaging over all other profileration rates. These visualisations are a way of understanding the full joint posterior distribution which is five-dimensional in full generality. Fig B. Posterior for proliferation rates under Ara-C treatment. The upper diagonal represents the marginal distributions for each proliferation rate when averaging over all other profileration rates. The plots below the diagonal show bivariate marginal distributions illustrating pairwise associations after averaging over all other profileration rates. These visualisations are a way of understanding the full joint posterior distribution which is five-dimensional in full generality. Fig C. Posterior for proliferation rates when recovering from Ara-C treatment. The upper diagonal represents the marginal distributions for each proliferation rate when averaging over all other profileration rates. The plots below the diagonal show bivariate marginal distributions illustrating pairwise associations after averaging over all other profileration rates. These visualisations are a way of understanding the full joint posterior distribution which is five-dimensional in full generality. Fig D. Typical section obtained during the experimental procedures described in the main manuscript. These are also described more fully in the companion paper [18].



OJM would like to thank the members of the Wolfson Centre of Mathematical Biology for useful discussion and feedback, in particular Ruth Baker for her helpful comments on stochastic modelling.

Funding Statement

This work was funded by The Biotechnology and Biological Sciences Research Council (, project numbers BB/K018256/1 (PKM, HMB, CP, OJM, AP), BB/K017578/1 (PKM, HMB, OJM), BB/K017144/1 (SRC, AJMW, AP), BB/J004529/1 (SRC, AJMW, CP, AP) and the The Engineering and Physical Sciences Research Council ( project number EP/I017909/1 (PKM, AGF). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability

Data Availability

All manuscript, code and data files can be found at the following github page:


1. Wright NA, Alison M (1984) The biology of epithelial cell populations. Oxford University Press, USA
2. Radtke F, Clevers H (2005) Self-renewal and cancer of the gut: Two sides of a coin. Science 307:1904–1909 doi: 10.1126/science.1104815 [PubMed]
3. Reuss L (2010) Epithelial transport Comprehensive physiology
4. van der Flier LG, Clevers H (2009) Stem cells, self-renewal, and differentiation in the intestinal epithelium. Annu Rev Physiol 71:241–260 doi: 10.1146/annurev.physiol.010908.163145 [PubMed]
5. Turner JR (2009) Intestinal mucosal barrier function in health and disease. Nat Rev Immunol 9:799–809 doi: 10.1038/nri2653 [PubMed]
6. Marchiando AM, Graham WV, Turner JR (2010) Epithelial barriers in homeostasis and disease. Annu Rev Pathol 5:119–144 doi: 10.1146/annurev.pathol.4.110807.092135 [PubMed]
7. Barker N (2014) Adult intestinal stem cells: Critical drivers of epithelial homeostasis and regeneration. Nat Rev Mol Cell Biol 15:19–33 doi: 10.1038/nrm3721 [PubMed]
8. Kaur P, Potten CS (1986) Cell migration velocities in the crypts of the small intestine after cytotoxic insult are not dependent on mitotic activity. Cell Tissue Kinet 19:601–610 [PubMed]
9. Kaur P, Potten CS (1986) Circadian variation in migration velocity in small intestinal epithelium. Cell Tissue Kinet 19:591–599 [PubMed]
10. Tsubouchi S (1983) Theoretical implications for cell migration through the crypt and the villus of labelling studies conducted at each position within the crypt. Cell Tissue Kinet 16:441–456 [PubMed]
11. Dunn S-J, Näthke IS, Osborne JM (2013) Computational models reveal a passive mechanism for cell migration in the crypt. PLoS One 8:e80516 doi: 10.1371/journal.pone.0080516 [PMC free article] [PubMed]
12. Meineke FA, Potten CS, Loeffler M (2001) Cell migration and organization in the intestinal crypt using a lattice-free model. Cell Prolif 34:253–266 doi: 10.1046/j.0960-7722.2001.00216.x [PubMed]
13. Loeffler M, Stein R, Wichmann HE, Potten CS, Kaur P, Chwalinski S (1986) Intestinal cell proliferation. I. A comprehensive model of steady-state proliferation in the crypt. Cell Tissue Kinet 19:627–645 [PubMed]
14. Loeffler M, Potten CS, Paulus U, Glatzer J, Chwalinski S (1988) Intestinal crypt proliferation. II. Computer modelling of mitotic index data provides further evidence for lateral and vertical cell migration in the absence of mitotic activity. Cell Tissue Kinet 21:247–258 [PubMed]
15. Vyshemirsky V, Girolami MA (2008) Bayesian ranking of biochemical system models. Bioinformatics 24:833–839 doi: 10.1093/bioinformatics/btn475 [PubMed]
16. Kirk P, Thorne T, Stumpf MP (2013) Model selection in systems and synthetic biology. Current opinion in biotechnology 24:767–774 doi: 10.1016/j.copbio.2013.03.012 [PubMed]
17. Weavers H, Liepe J, Sim A, Wood W, Martin P, Stumpf MP (2016) Systems analysis of the dynamic inflammatory response to tissue damage reveals spatiotemporal properties of the wound attractant gradient. Current Biology 26:1975–1989 doi: 10.1016/j.cub.2016.06.012 [PMC free article] [PubMed]
18. Parker A, Maclaren OJ, Fletcher AG, Muraro D, Kreuzaler PA, Byrne HM, Maini PK, Watson AJM, Pin C (2016) Cell proliferation within small intestinal crypts is the principal driving force for cell migration on villi. The FASEB Journal; published ahead of print October 20, 2016. [PMC free article] [PubMed]
19. Bernardo JM, Smith AFM (2009) Bayesian theory. John Wiley & Sons
20. Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB (2013) Bayesian data analysis, third edition Taylor & Francis
21. Tarantola A (2005) Inverse problem theory and methods for model parameter estimation. SIAM
22. Berliner LM (1996) Hierarchical Bayesian time series models In: Maximum entropy and Bayesian methods. Springer; Netherlands, pp 15–22
23. Cressie N, Wikle CK (2011) Statistics for spatio-temporal data. John Wiley & Sons
24. Wikle CK (2015) Modern perspectives on statistics for spatio-temporal data. WIREs Comput Stat 7:86–98 doi: 10.1002/wics.1341
25. Dawid AP (2002) Influence diagrams for causal modelling and inference. Int Stat Rev 70:161–189 doi: 10.2307/1403901
26. Dawid PA (2010) Seeing and doing: The Pearlian synthesis In: Dechter Rina, Geffner Hector, and Halpern Joseph Y. (ed) Heuristics, probability and causality: A tribute to judea pearl. College Publications; London, pp 309–325
27. Dawid PA (2010) Beware of the DAG! NIPS Causality: Objectives and Assessment 6:59–86
28. Box GEP (1976) Science and statistics. J Am Stat Assoc 71:791–799 doi: 10.1080/01621459.1976.10480949
29. Box GEP (1980) Sampling and Bayes’ inference in scientific modelling and robustness. J R Stat Soc Ser A 143:383–430 doi: 10.2307/2982063
30. Evans M (2015) Measuring statistical evidence using relative belief. CRC Press [PMC free article] [PubMed]
31. Gelman A, Shalizi CR (2013) Philosophy and the practice of Bayesian statistics. Br J Math Stat Psychol 66:8–38 doi: 10.1111/j.2044-8317.2011.02037.x [PMC free article] [PubMed]
32. Kozar S, Morrissey E, Nicholson AM, van der Heijden M, Zecchini HI, Kemp R, Tavaré S, Vermeulen L, Winton DJ (2013) Continuous clonal labeling reveals small numbers of functional stem cells in intestinal crypts and adenomas. Cell Stem Cell 13:626–633 doi: 10.1016/j.stem.2013.08.001 [PubMed]
33. Vermeulen L, Morrissey E, van der Heijden M, Nicholson AM, Sottoriva A, Buczacki S, Kemp R, Tavaré S, Winton DJ (2013) Defining stem cell dynamics in models of intestinal tumor initiation. Science 342:995–998 doi: 10.1126/science.1243148 [PubMed]
34. Lopez-Garcia C, Klein AM, Simons BD, Winton DJ (2010) Intestinal stem cell replacement follows a pattern of neutral drift. Science 330:822–825 doi: 10.1126/science.1196236 [PubMed]
35. Meineke FA, Potten CS, Loeffler M (2001) Cell migration and organization in the intestinal crypt using a lattice-free model. Cell Prolif 34:253–266 doi: 10.1046/j.0960-7722.2001.00216.x [PubMed]
36. Potten CS, Roberts SA, Chwalinski S, Loeffler M, Paulus U (1988) Scoring mitotic activity in longitudinal sections of crypts of the small intestine. Cell Tissue Kinet 21:231–246 [PubMed]
37. Pearl J (2009) Causality. Cambridge University Press
38. Pearl J (2009) Causal inference in statistics: An overview. Stat Surv 3:96–146 doi: 10.1214/09-SS057
39. Dawid PA (1979) Conditional independence in statistical theory. J R Stat Soc Series B Stat Methodol 41:1–31
40. Woodward J (2003) Making things happen: A theory of causal explanation. Oxford University Press
41. Woodward J (1997) Explanation, invariance, and intervention. Philos Sci 64:S26–S41 doi: 10.1086/392584
42. Spirtes P, Glymour CN, Scheines R (2000) Causation, prediction, and search. MIT Press
43. Van Kampen NG (1992) Stochastic processes in physics and chemistry. Elsevier
44. Maclaren OJ, Byrne HM, Fletcher AG, Maini PK (2015) Models, measurement and inference in epithelial tissue dynamics. Math Biosci Eng 12:1321–1340 doi: 10.3934/mbe.2015.12.1321 [PubMed]
45. Wilkinson DJ (2011) Stochastic modelling for systems biology. CRC Press
46. Baker RE, Yates CA, Erban R (2010) From microscopic to macroscopic descriptions of cell migration on growing domains. Bull Math Biol 72:719–762 doi: 10.1007/s11538-009-9467-x [PubMed]
47. Hywood JD, Hackett-Jones EJ, Landman KA (2013) Modeling biological tissue growth: Discrete to continuum representations. Phys Rev E Stat Nonlin Soft Matter Phys 88:032704 doi: 10.1103/PhysRevE.88.032704 [PubMed]
48. Fozard JA, Byrne HM, Jensen OE, King JR (2010) Continuum approximations of individual-based models for epithelial monolayers. Math Med Biol 27:39–74 doi: 10.1093/imammb/dqp015 [PubMed]
49. Murray PJ, Edwards CM, Tindall MJ, Maini PK (2009) From a discrete to a continuum model of cell dynamics in one dimension. Phys Rev E Stat Nonlin Soft Matter Phys 80:031912 doi: 10.1103/PhysRevE.80.031912 [PubMed]
50. Rasmussen C, Williams C (2005) Gaussian processes for machine learning. MIT
51. Johnston MD, Edwards CM, Bodmer WF, Maini PK, Chapman SJ (2007) Examples of mathematical modeling: Tales from the crypt. Cell Cycle 6:2106–2112 doi: 10.4161/cc.6.17.4649 [PMC free article] [PubMed]
52. Carulli AJ, Samuelson LC, Schnell S (2014) Unraveling intestinal stem cell behavior with models of crypt dynamics. Integr Biol 6:243–257 doi: 10.1039/c3ib40163d [PMC free article] [PubMed]
53. Robert C, Casella G (2013) Monte carlo statistical methods. Springer Science & Business Media
54. Ketcheson DI, Mandli K, Ahmadia AJ, Alghamdi A, de Luna MQ, Parsani M, Knepley MG, Emmett M (2012) PyClaw: Accessible, extensible, scalable tools for wave propagation problems. SIAM J Sci Comput 34:C210–C231 doi: 10.1137/110856976
55. Mandli KT, Ketcheson DI, others (2014) PyClaw software.
56. Clawpack Development Team (2014) Clawpack software.
57. Guyer JE, Wheeler D, Warren JA (2009) FiPy: Partial differential equations with python. Comput. Sci. Eng. doi: 10.1109/MCSE.2009.52
58. Gelman A (2004) Exploratory data analysis for complex models. J Comput Graph Stat 13:755–779 doi: 10.1198/106186004X11435
59. LeVeque RJ (2002) Finite volume methods for hyperbolic problems. Cambridge University Press
60. Aguilar O, Allmaras M, Bangerth W, Tenorio L (2015) Statistics of parameter estimates: A concrete example. SIAM Rev 57:131–149 doi: 10.1137/130929230
61. Kershaw SK, Byrne HM, Gavaghan DJ, Osborne JM (2013) Colorectal cancer through simulation and experiment. IET Syst Biol 7:57–73 doi: 10.1049/iet-syb.2012.0019 [PubMed]
62. De Matteis G, Graudenzi A, Antoniotti M (2013) A review of spatial computational models for multi-cellular systems, with regard to intestinal crypts and colorectal cancer development. J Math Biol 66:1409–1462 doi: 10.1007/s00285-012-0539-4 [PubMed]
63. Fletcher AG, Murray PJ, Maini PK (2015) Multiscale modelling of intestinal crypt organization and carcinogenesis. Math Models Methods Appl Sci 25:2563–2585 doi: 10.1142/S0218202515400187
64. Berliner LM (2003) Physical-statistical modeling in geophysics. J. Geophys. Res. 108: doi: 10.1029/2002JD002865
65. Wikle CK (2003) Hierarchical models in environmental science. Int Stat Rev 71:181–199 doi: 10.1111/j.1751-5823.2003.tb00192.x
66. Cressie N, Calder CA, Clark JS, Ver Hoef JM, Wikle CK (2009) Accounting for uncertainty in ecological analysis: The strengths and limitations of hierarchical statistical modeling. Ecol Appl 19:553–570 doi: 10.1890/07-0744.1 [PubMed]
67. Ogle K (2009) Hierarchical Bayesian statistics: Merging experimental and modeling approaches in ecology. Ecol Appl 19:577–581 doi: 10.1890/08-0560.1 [PubMed]
68. Blei DM (2014) Build, compute, critique, repeat: Data analysis with latent variable models. Annual Review of Statistics and Its Application 1:203–232 doi: 10.1146/annurev-statistics-022513-115657
69. Pawitan Y (2001) In all likelihood: Statistical modelling and inference using likelihood. OUP Oxford
70. Pawitan Y, Lee Y (2016) Wallet game: Probability, likelihood and extended likelihood. Am Stat 0:1–7
71. Lee Y, Nelder JA, Pawitan Y (2006) Generalized linear models with random effects: Unified analysis via h-likelihood. CRC Press
72. Mosegaard K, Tarantola A (2002) Probabilistic approach to inverse problems. International Geophysics Series 81:237–268 doi: 10.1016/S0074-6142(02)80219-4
73. Wood SN (2010) Statistical inference for noisy nonlinear ecological dynamic systems. Nature 466:1102–1104 doi: 10.1038/nature09319 [PubMed]
74. Marin J-M, Pudlo P, Robert CP, Ryder RJ (2012) Approximate Bayesian computational methods. Stat Comput 22:1167–1180 doi: 10.1007/s11222-011-9288-2
75. Wilkinson RD (2013) Approximate Bayesian computation (ABC) gives exact results under the assumption of model error. Stat Appl Genet Mol Biol 12:129–141 doi: 10.1515/sagmb-2013-0010 [PubMed]
76. Ratmann O, Andrieu C, Wiuf C, Richardson S (2009) Model criticism based on likelihood-free inference, with an application to protein network evolution. Proc Natl Acad Sci U S A 106:10576–10581 doi: 10.1073/pnas.0807882106 [PubMed]
77. Davies PL (2014) Data analysis and approximate models: Model choice, Location-Scale, analysis of variance, nonparametric regression and image analysis. CRC Press

Articles from PLoS Computational Biology are provided here courtesy of Public Library of Science