Modelling cell turnover in a complex tissue during development

Journal of Theoretical Biology 338 (2013) 66–79

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Modelling cell turnover in a complex tissue during development J. Lefevre, D.J. Marshall, A.N. Combes, A.L. Ju, M.H. Little, N.A. Hamilton n Institute for Molecular Bioscience, The University of Queensland, St. Lucia, Brisbane, QLD 4072, Australia

H I G H L I G H T S







We develop a model of cell cycle and proliferation in an organ sub-compartment during development. The model combined with experimental imaging predicts factors driving compartment dynamics. The method is demonstrated on experimental data from the cap mesenchyme of the developing kidney. Significant changes in cell cycle length and cell differentiation are detected by the model.

art ic l e i nf o

a b s t r a c t

Article history: Received 7 December 2012 Received in revised form 26 August 2013 Accepted 28 August 2013 Available online 7 September 2013

The growth of organs results from proliferation within distinct cellular compartments. Organ development also involves transitions between cell types and variations in cell cycle duration as development progresses, and is regulated by a balance between entry into the compartment, proliferation of cells within the compartment, acquisition of quiescence and exit from that cell state via differentiation or death. While it is important to understand how environmental or genetic alterations can perturb such development, most approaches employed to date are descriptive rather than quantitative. This is because the identification and quantification of such parameters, while tractable in vitro, is challenging in the context of a complex tissue in vivo. Here we present a new framework for determining cell turnover in developing organs in vivo that combines cumulative cell-labelling and quantification of distinct cell-cycle phases without assuming homogeneity of behaviour within that compartment. A mathematical model is given that allows the calculation of cell cycle length in the context of a specific biological example and assesses the uncertainty of this calculation due to incomplete knowledge of cell cycle dynamics. This includes the development of a two population model to quantify possible heterogeneity of cell cycle length within a compartment and estimate the aggregate proliferation rate. These models are demonstrated on data collected from a progenitor cell compartment within the developing mouse kidney, the cap mesenchyme. This tissue was labelled by cumulative infusion, volumetrically quantified across time, and temporally analysed for the proportion of cells undergoing proliferation. By combining the cell cycle length predicted by the model with measurements of total cell population and mitotic rate, this approach facilitates the quantification of exit from this compartment without the need for a direct marker of that event. As a method specifically designed with assumptions appropriate to developing organs we believe this approach will be applicable to a range of developmental systems, facilitating estimations of cell cycle length and compartment behaviour that extend beyond simple comparisons of mitotic rates between normal and perturbed states. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Cell cycle Cell differentiation Modelling Kidney Development

1. Introduction Being able to determine the magnitude and dynamics of a cell population over time is invaluable for analyzing tissue development, homeostasis, and the origins of disease. Most cell populations do not occupy a steady state but instead are in flux, notably so in developing tissues. For any given cell population within a tissue, the size of the population will depend upon the rate of

n

Corresponding author. Tel.: þ 61 7 3346 2033; fax: þ61 7 3346 2101. E-mail address: [email protected] (N.A. Hamilton).

0022-5193/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jtbi.2013.08.033

proliferation, entry into or exit of cells from that population where entry can include inward migration or other cells adopting that particular fate and exit can include cell death or differentiation into a distinct cell type (Fig. 1D). Proliferation involves the passage of a cell through the phases of the cell cycle to generate two daughter cells (Fig. 1E). While markers are available to identify what phase of the cell cycle a given cell is passing through, the length of each phase can vary and cells within a given population may be actively passing through the cell cycle or resting in G0 (quiescence) (Fig. 1E) (Cooper, 2000). For those cells actively proliferating, a variety of regulators of gap phases G1 and G2 can affect the rate at which a given cell passes through the cell cycle

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

67

Fig. 1. (A) Optical section of a whole mouse embryonic kidney at 15.5 dpc. Cap mesenchyme cells are highlighted in green with an antibody to Six2, cells undergoing mitosis are in red (antibody to Phh3), and the ureteric tree is labelled in white (antibody to Calbindin), scale bar 200 microns. (B) Zoom of the region outlined in A. The cap mesenchyme (CM) compartment is outlined, as is the early nephron compartment (N). (C) A cartoon of the CM compartment (CM) and exit of cells from this compartment to the early nephron (N) via differentiation (arrow). (D) Turnover within a specific cellular compartment. Diagram indicates the parameters determining change in compartment size and the end result of different types of proliferation. All cells within the box (cell compartment of interest) are identified as remaining in this state based upon expression of a specific marker (Six2 in the case of the cap mesenchyme of the kidney). Note that the entry of a cell in the compartment into a quiescent state (G0) is not regarded as reducing the size of that compartment as these cells are assumed to remain marker positive. Exit from the compartment is viewed as loss of that cell state and may occur via cell death or differentiation into another cell state (as defined by loss of the compartment-specific marker). (E) Diagram of the cell cycle indicating markers for S-phase (S) and mitosis (M) and showing the chronology of the gap phases G1, G2 and Quiescence (G0). Exit from the cell cycle occurs either by entry into G0, which is potentially reversible, or by cell death, and should not be confused with exit from the compartment of interest via differentiation into another state, shown in panel D. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

and cell division within a population is usually asynchronous (Cooper, 2000). Hence, the average proportion of cells within any one phase of the cell cycle can only yield relative proliferation comparisons and not absolute measures of cell cycle length. The passage of any cell through the cell cycle can also yield different outcomes. Many developing organs contain progenitor populations responsible for generating more differentiated cell types. These populations are often self-renewing, however this can include both asymmetric self-renewal (parent cell replaces itself and generates one distinct progeny) and symmetric self-renewal (parent cell replaces itself with two other cells of the same type) (Fig. 1D). Even symmetric cell division can result in the formation of two differentiated cells as opposed to doubling the number of the initial cell type (Fig. 1D). Without a marker for the differentiated state and/or a lineage tracing method for showing what cell gave rise to what progeny, it may not be possible to document exit via differentiation. Finally, cell death can result in a decline in population that may artificially appear as a reduced level of proliferation (Fig. 1D). As a result, modelling the turnover of a cellular population, even in vitro, is complicated. This is compounded when the population of interest

is within a complex tissue or developing organ. Note that in the current paper we use exit to refer to exit from the compartment, rather than exit from the cell cycle via quiescence or death (Fig. 1). In order to fully understand defects in development or disease, we require a capacity to morphometrically and quantitatively model the rate of cell turnover within complex tissues. Early work (Steel and Bensted, 1965; Steel, 1967, 1977) in this area dealt with tumour cell dynamics. A more recent paper is this area (Bertuzzi et al., 1997) provides a detailed analysis of cell cycle dynamics in tissue affected by non-uniform loss, but does not consider cumulative S-phase marker experiments. A large body of more recent work (see (Nowakowski et al., 1989; Takahashi et al., 1993, 1995; Caviness and Takahashi, 1995; Cai et al., 1997; Kornack and Rakic, 1998; Calegari et al., 2005; Charvet and Striedter, 2008) and many other papers citing (Nowakowski et al., 1989)) addresses cell cycle length in the developing neocortex and other brain tissue, using cumulative cell cycle phase labelling strategies as well as multiple labelling techniques, whilst in (Alexiades and Cepko, 1996), quiescent cells are also considered and multiple markers used inside the developing rat retina. These studies involve isolated populations

68

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

and often ignore the effect of entry and exit, important considerations particularly during development and repair. However the most critical issue is that to our knowledge these methods all use the linear model given in (Nowakowski et al., 1989) for the interpretation of cumulative S-phase marker experiments, which assumes a steady state homogenous non-quiescent population in which the total volume or number of cells in a tissue is assumed to be constant with proliferation simply replacing those cells that have exited though processes such as apoptosis. In a developing organ such as the kidney, there is often rapid proliferation with cell numbers increasing by orders of magnitude in a matter of days (see Section 3.1 below). Other studies of cell populations over time involve mathematical models with many parameters, as in (Baker et al., 1998), not all of which may be accessible. Fitting such models to experimental data can be problematic in that solutions may not be unique or are sensitive to initial conditions and care must be taken when fitting them to data. Multiple sources of information are often necessary to reduce the number of free parameters. When studying development in the mouse kidney it became apparent that there was a need to develop an understanding of and explicitly model the full range of cell behaviours inside selfrenewing compartments over time. Our goal was to understand what drives change in population size during organogenesis and to be able to mathematically model the outcome for the organ, in particular when direct markers for certain behaviours are absent. Here we develop a cell cycle and proliferation model appropriate to the timescales that occur in development in an identifiable compartment of an organ. The compartment represents a progenitor population within a developing organ which is growing in size as well as generating new cell types via differentiation. The method brings together mathematical modelling and quantified experimental data to provide a flexible framework to study cell turnover. Experimentally, a cumulative marker of cell division is central to the approach to generate data on proliferation within a cellular compartment. In this study we use the thymidine analog EdU. A mathematical model of the expected proportion of marked to unmarked cells over time is developed and from parameter fitting to the experimental data, cell cycle length and S-phase duration are determined. Care is taken to specify assumptions clearly and not introduce parameters unnecessarily. Key elements of the model are that homeostasis is not assumed and cells are asynchronous in their cell cycles. The strategy is then to combine the cell cycle information so determined with timed, static measurements of compartment volume and proportion of cells in mitosis. As we show, this enables the extrapolation to earlier and later time points and provides information about changes in proliferation or exit from that population across development without the requirement for explicit markers of differentiation. These combined modelling and experimental methods are demonstrated using the nephron progenitor population inside the developing kidney. This population, also referred to as the cap mesenchyme, is defined by the presence of Cited1 and Six2 proteins and is a self-renewing population of progenitor cells that differentiate into the epithelial nephrons of the kidney through a mesenchyme to epithelial transition(Georgas et al., 2009; Boyle et al., 2008; Kobayashi et al., 2008; Hendry et al., 2011). Previous studies have shown using lineage tracing that all differentiated cells within the nephrons, aside from the collecting duct, arise from cap mesenchyme cells and that the starting size of this population is determined at approximately 11.5 days of development (days post coitum; dpc) in mice (Kobayashi et al., 2008). This defines the population as one in which there is no entry, leaving only exit (death and differentiation) and proliferation to regulate the size of the population over time. See, for example, (Hendry et al., 2011) for a detailed description of the biological context and key developmental role of the cap mesenchyme. Using the

methods developed in this paper, we define key parameters of cap mesenchyme morphogenesis such as cell cycle length and S-phase duration, and we have quantified the level of cell exit via differentiation from the compartment during subsequent development.

2. Methods The aim of this work is to understand and model cell turnover in a marked compartment by quantifying each of the cell behaviours that contribute to change in the population size over time. In the following, we begin by describing the biological assumptions about the proliferating populations under study (Section 2.1). Primarily we will be concerned with the cap mesenchyme in the developing kidney, though the assumptions made should be applicable to a wide range of developing organs and cell populations. Subsequently, in Section 2.2, we develop mathematical models based on these biological assumptions to enable quantitative interpretation of experimental data and prediction of quantities such as cell cycle length. 2.1. Biological assumptions and experiments The cell population behaviours to be considered are entry, proliferation, quiescence, and exit from the compartment via death or differentiation as shown in Fig. 1D. 2.1.1. Proliferation, exit and entry Proliferation occurs when cell spass through the cell cycle and undergo cell division. Quiescence describes cells that have exited the cell cycle (but not necessarily the compartment) and hence are not part of the proliferating population. A cell may leave the compartment of interest via either cell death (generally apoptosis) or differentiation into a different cell type, collectively referred to as exit. Entry refers to the transition of a cell of another type into the compartment. Here, entry into G0 (quiescence) is not described as exit as such cells would remain positive for the compartment-specific marker used to quantify the size of that compartment. A representation of the cell cycle is shown in Fig. 1E. When a cell divides, the progeny can either become dormant (quiescent), die via apoptosis, or re-enter the cell cycle and divide again at mitosis. Surviving progeny may either remain within the compartment or differentiate. Division may be symmetric, with the progeny of a cell sharing the same fate, or asymmetric. The cap mesenchyme is known to contain a self-renewing cell population as well as exhibiting exit via differentiation into renal vesicle (RV) cells, an event triggered by a specific growth factor and accompanied by a loss of compartment-specific markers (Kobayashi et al., 2008). Cell death and quiescence are not known to occur, but must be considered as possibilities. In the case of the developing cap mesenchyme, lineage studies (Kobayashi et al., 2008) do indicate that no entry occurs after initial formation of the compartment (although in general cases such behaviour can be considered). Thus entry is the only dynamic excluded a priory. The mechanism of exit from the CM is not known with certainty, so in addition to the model of exit occurring immediately after mitosis, we also consider the possibility that exit may occur without regard to the stage of the cell cycle, for example in response to a positional signal. Upon differentiation into RV, the CM marker (Six2) is lost. Proliferation rates and cell counts in the developing nephron are unknown, and so do not provide any information on exit via differentiation. There is also no direct marker of differentiation available. Therefore exit can only be inferred from quantitation of the CM, as shall be shown later by combining experimental results with cell cycle modelling.

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

2.1.2. Cell cycle and synchronisation Cell cycle length is known to vary during development (Cooper, 2000) due to differences in requirements for individual cell growth and morphological complexity and the requirements for gap phases for repair. We assume that any change is sufficiently slow that the cell cycle length can be regarded as constant over an experimental timecourse of less than one cell cycle. In the first instance we also assume that at a particular point of time each phase of the cell cycle is of fixed length across all cells in the compartment. In order to allow for the possibility of heterogeneity of proliferation while minimising the number of additional parameters, we also later consider a model in which the compartment contains two subpopulations, differing only in the length of the first gap phase. For typical lengths of the phases see (Cooper, 2000). While all cells of a particular state progress through division at the same rate, they usually do not do this in synchrony, that is, simultaneous mitosis, division and so forth. Hence, we assume asynchrony in the observed population of cells progressing through the cell cycle, meaning that the fraction of cells in any particular phase of the cell cycle is independent of time. However we also consider the possibility of an element of random synchrony in individual samples. This will be discussed further in Section 2.2.1. 2.1.3. Compartment and cell cycle markers In order to detect and quantify cellular state, a number of experimental markers are used and imaged. Cells were identified with a compartment-specific marker and co-localised with markers of particular cell phases. For information on the specific markers used for the experiments, see Appendix 5.1. In brief, the cap mesenchyme population was identified via fluorescent labelling with an antibody to the Six2 protein, a transcription factor specific to the cap mesenchyme,

69

as described in Appendix 5.1. An example of such imaging can be seen in Fig. 2C. A cumulative marker, EdU, was used to mark cells as they entered into S-phase, with marked cells retaining the marker over time due to stable DNA intercalation during S-phase. Fig. 2A shows Six2 cells highlighted in green, while EdU positive cells are shown in the red channel. Thus, yellow regions delineate double labelled cells, cap mesenchyme cells that are in or have passed through S-phase. A static mitosis marker, pHH3, and a marker for death via apoptosis, caspase, were also used to detect cells in these states. Fig. 2D shows Six2 positive (cap mesenchyme) cells in green with apoptotic caspase positive cells in red. It is worth noting that few, if any, of the apoptotic cells appear to be part of the cap mesenchyme. It will be assumed that apoptosis is the only form of cell death in the compartment and it is effectively absent. In Fig. 2B, Six2 positive cells (cap mesenchyme) are again shown in green, with pHH3 positive (mitotic) cells shown in the red channel. Again, coincident red and green (yellow) shows those cap mesenchyme cells that were undergoing mitosis. In the following section we will consider how the combination of these experimental data types, together with appropriate models, may be utilised to extract information on quiescence and the length of the cell cycle length over time in the cap mesenchyme. 2.2. Modelling division and quiescence with a cumulative S-phase marker In this section we consider the appropriate approach to use in the modelling and interpretation of cumulative S-phase marker data in the biological context of the CM as outlined above, and also in the context of the extensive prior work in this area. In particular

Fig. 2. Representative slices from image data. (A) Confocal image of Six2 þ (green) and EdU þ (red) cells within the cap mesenchyme. Colocalised signal is yellow and marks double positive cells (cap mesenchyme cells labelled for EdU). Arrows show some of the many double positive cells. (B) Confocal image of Six2 þ (green) and pHH3 þ (red) cells, with calbindin (white) marking ureteric tip epithelium. Arrows show examples of double positive Six2 þ /pHH3þ cells. (C) Slice through optical projection tomography (OPT) image of whole kidney, Six2 channel only. (D) Confocal image of Six2þ (green) and caspase þ (red) cells within the cap mesenchyme, revealing a paucity of cells undergoing apoptosis within this compartment. The only caspase-positive cells present (arrowheads) are found within renal vesicle structures and not in the Six2 þ (green) cap mesenchyme compartment. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

70

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

we develop an approach to account for some areas of uncertainty regarding exit and the possibility of cell cycle synchrony. 2.2.1. General model In the following, we assume the cell population of interest is distinguished by a marker (such as Six2) and that a cumulative S-phase marker (EdU) is introduced at a time t ¼ 0, where t is in hours, with booster injections to ensure that any cell that enters S-phase during the period of observation will take up and exhibit the S-phase marker. Observation is on a proportion of the compartment only, so the experimental measurement is of the proportion of EdU þ cells rather than absolute numbers. For more information on experimental methods, see Appendix 5.1. To study a cumulative S-phase marker over time, a model for cell proliferation needs to be considered. Models of cell proliferation over time have been studied extensively in the literature (see (Baker et al., 1998) for a good survey). An important consideration is the distribution of cell cycles at a particular point in time. We initially assume that the length of the cell cycle and the phases within it are constant within the CM compartment at a given developmental stage, with no quiescence. Let T G1 , T S , T G2 and T m be the durations of G1, S-phase, G2 and mitosis respectively, and T c ¼ T G1 þT S þ T G2 þ T m be the overall cell cycle length. The difference between the G2 and mitosis phases is irrelevant to the EdU experiment, so we let T G2 þ m ¼ T G2 þ T m . We define the age of an actively cycling cell as the time since the previous mitosis event and entry into G1 (note that the current cell cycle phase of a given cell is a simple function of age under these assumptions). Let nða; tÞ be the cell density with respect to age a at time t; that is, Ra at time t the number of cells aged between a1 and a2 is a12 nða; tÞda. RT The total number of cells at time t is given by NðtÞ ¼ 0 c nða; tÞda. All cells that are or were in S-phase at any point of the EdU timecourse are assumed to be EdU þ , as will the progeny of such cells. Although mitosis of an EdU cell will result in weaker EdU expression in the daughter cells, regular EdU delivery will mean that there will be at most one such dilution event, which is assumed to leave detectable EdU expression. Therefore the number of EdUþ cells is 8 R T þT þt G1 S > nða; tÞda; 0 rt o T G2 þ m ; > < T G1 R TC R tT EðtÞ ¼ nða; tÞda þ 0 G2 þ m nða; tÞda; T G2 þ m r t o T c T s : T G1 > > : NðtÞ; t Z T T : c

s

Since we are limited to counting EdU þ and EdU  cells within a sample of the compartment, the measured quantity is the EdU proportion RðtÞ ¼ EðtÞ=NðtÞ. The standard assumption of asynchrony, that the proportion of cells at each stage of the cell cycle is constant over time, is equivalent to the condition that nða; tÞ is multiplicatively separable (Fig. 3). The simplest asynchronous model is that nða; tÞ is constant with respect to a; that is, that all ages are equally likely. This rectangular distribution

gives RðtÞ ¼ ðT s þ tÞ=T c for 0 r t r T c T s ; and RðtÞ ¼ 1 thereafter. However, this linear model can only be justified under a steadystate assumption, in which on average each mitotic event results in exactly one of the daughter cells re-entering the population, with the other exiting through death or differentiation. In the absence of exit (and entry), it is necessary that ðd=dtÞnða; tÞ ¼ ðd=daÞnða; tÞ, due to cells progressing through the cycle, and nð0; tÞ ¼ 2  nðT c ; tÞ, due to mitosis. The unique asynchronous distribution satisfying these constraints is the exponential model nða; tÞ ¼ k  2ðtaÞ=T c , for some constant k (for instance see (Bonhoeffer et al., 2000)). The linear model has also been used where the steady state assumption is considered reasonable, or for simplicity (see, for instance, (Nowakowski et al., 1989; Kornack and Rakic, 1998; Calegari et al., 2005); Charvet and Striedter, 2008)). A complicating factor in the CM is exit from the compartment via differentiation. Although the CM compartment grows rapidly, it is known that some degree of exit through differentiation occurs; quantification of this is a principal goal of this paper. The situation is analogous to that studied in (Bertuzzi et al., 1997) in the context of cell death in human tumours. If exit occurs at any point in the cell cycle without bias with respect to age, then the proportion of cells at each age will be unaffected. Although the above formula for nða; tÞ will not strictly hold, the value of RðtÞ is unchanged and hence such exit may be ignored. However, it is considered likely that differentiation occurs after mitosis, perhaps due to asymmetrical division of progenitor cells. Suppose that the proportion of daughter cells remaining within the compartment is θ. Then the relationship between the density function before and after mitosis is modified to nð0; tÞ ¼ 2θ  nðT c ; tÞ, implying the asynchronous distribution is nða; tÞ ¼ k  r ðtaÞ=T c , where r ¼ 2θ. This is a generalisation of the standard exponential distribution, and reduces to the linear model if r ¼ 1. Performing the appropriate integrations for EðtÞ and NðtÞ and dividing gives the formula 8 ðT G2 þ m þ T s Þ=T c r ðT G2 þ m tÞ=T c Þ=ðr1Þ; 0 r t o T G2 þ m ; > < ðr RðtÞ ¼ 1ðr ðT c þ T G2 þ m tÞ=T c r ðT G2 þ m þ T s Þ=T c Þ=ðr1Þ; T G2 þ m r t o T c T s ; > : 1; t Z T c T s ;

for any positive r except r ¼ 1, in which case the linear model holds. Provided that the cell population is not shrinking, then r lies between 2 (if exit is negligible or occurs without age bias) and 1 (if exit occurs at mitosis and approaches steady-state levels). Hence the true value for RðtÞ is unknown but intermediate between the standard exponential model (above formula with r ¼ 2Þ and the linear model. Another factor to consider is that the assumption of asynchrony may not completely hold in practice, particularly in the case of an embryonic organ undergoing rapid proliferation, with cell behaviour evolving on a timescale comparable to the cell cycle length. Any synchrony is potentially complex, however we can estimate the maximum impact on the EdU results by considering the extreme case of complete synchronisation, in which all cells in the compartment have the same age. That is, for given t, nða; tÞ consists of a delta function at some age a0 . Now in practice, each EdU data point is

Fig. 3. For the non-quiescent population, the proportion of EdU þ cells is T S =T C at t¼ 0 if the distribution is uniform, however exponential growth leads to a non-uniform age distribution and more complex relationship, while synchronous cell cycles lead to variability in the proportion of cells in S-phase. The saturation time is T C T S provided that the cell cycle length is homogenous. Boxes represent S-phase and horizontal lines represent other elements of the cell cycle.

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

obtained from a separate embryo and the timing of developmental age is only precise to within a few hours, so a0 may be safely considered to be a uniform random variable on the interval ½0; T c , with separate and independent values for each sample. This implies high variability in measured EdU incorporation proportions, with individual samples exhibiting either no EdU incorporation or complete saturation of the cycling population, in this hypothetical extreme case. But the probability of a positive result (all cells EdU þ ) will be equal to ðT s þ tÞ=T c , so averaged results will recapitulate the linear model. Due to all the above considerations, the method used is to fit the standard exponential model (r ¼ 2 above) to obtain our primary result, but also to fit the linear model as a sensitivity test, giving an upper bound on the potential impact of synchrony and exit on the fitted parameters. In Section 3.1 this approach will be applied to experimental data to estimate T c and T s in the CM of a developing mouse embryo kidney. Note that T G2 þ m is a parameter in the exponential model (the rate of change in EdU incorporation jumps upwards at this time-point as existing EdUþ cells start to divide, adding to the increase caused by cells entering S-phase), but not in the linear model. 2.2.2. Quiescence The general model outlined above assumes that all cells in the defined compartment are actively engaged in the cell cycle (that is the growth fraction, the proportion of non-quiescent cells, is equal to 1). Allowing for quiescence introduces some complexities to the analysis, as cells (which may be EdU þ ) must exit from the cycling population into quiescence if the quiescent population is not to dwindle into insignificance. These issues have been treated in the literature (for example (Alexiades and Cepko, 1996)), but we take the simpler approach of attempting to continue the EdU time series to saturation, in order to demonstrate the absence of significant quiescence (at least at the start of the timecourse), and rendering such considerations superfluous. 2.2.3. Two population models The cap mesenchyme is known to exhibit heterogeneity with respect to marker gene expression (Mugford et al., 2009). Hence while all of the compartment is Six2 positive, this can be subsected on the expression of other genes. While this heterogeneity has been predicted to represent phenotypic heterogeneity, this has not been proven. However, it is reasonable to ask whether this may be reflected in heterogeneity of cell cycle length. Towards modelling this, we take the simplest plausible approach, assuming two distinct subpopulations which differ only in the length of the first gap phase. In the absence of entry to or exit from either subpopulation, the faster cycling subpopulation would become dominant over time. There is no indication that cellular heterogeneity in the cap mesenchyme is transient, so we do not favour a model implying that proliferative heterogeneity is transient. Instead we assume that the relative sizes of the subpopulations are approximately constant. One possible mechanism allowing this is a steady transition of cells from the faster subpopulation to the slower subpopulation. This introduces a complication also encountered when modelling quiescence (which can be regarded as a special case of the two population model), since cells entering a subpopulation will not in general resemble the existing members of that population in age or EdU incorporation. However, in the absence of direct evidence on this issue we considered the most plausible mechanism to be differential exit from the compartment, with little or no transition between the subpopulations. Disproportionate exit from the faster cycling subpopulation is assumed to maintain the relative subpopulation sizes. Under these assumptions, the EdU incorporation rates of the subpopulations each follow a

71

single population model, and the overall incorporation function is a weighted average of the incorporation functions of the two subpopulations, as used in Steel (1967). Let T C1 and T C2 be the cell cycle lengths of the two subpopulations, T C1 o T C2 , and let a be the proportion of cells in the slower cycling subpopulation. By assumption the subpopulations have the same length S-phase (T s ). Assuming the subpopulations follow the asynchronous exponential model with G2 and mitosis having negligible length, then 8 ðT c1 tÞ=T c1 2T s =T c1 Það2ðT c2 tÞ=T c2 2T s =T c2 Þ; 0 r t oT C1 T s ; > < 1ð1aÞð2 ðT c2 tÞ=T c2 T s =T c2 Rðt Þ ¼ 1að2 2 Þ; T C1 T s r t o T C2 T s ; > : 1; t Z T c2 T s :

In the case of a two population model, it is useful to give an effective overall cell cycle length which is comparable to the single population model cell cycle length; we use the weighted harmonic mean T C ¼ 1=ðð1a=T C1 Þ þ ða=T C2 ÞÞ. This is the cell cycle length for a homogenous population with the same rate of proliferation, or equivalently the time taken for the population to double in the absence of entry, exit or change in the rate of proliferation; this is Steel's potential doubling time (see for example (Steel, 1967; Bertuzzi et al., 1997)).

2.3. Mitosis and apoptosis As described in Section 5.1, experimental markers such as phosphohistone H3 may be used to mark those cells undergoing mitosis at a given time point. Expected fractions of cells marked as being in mitosis can then be determined much as for a S-phase marker. Let t be time from the start of the timecourse and suppose the population follows an asynchronous exponential distribution with no quiescent cells in the population. Given that T m {T c , the proportion of the cycling cell population that is in mitosis is closely approximated by T m nðT c ; tÞ=ðT c NðtÞÞ  lnð2ÞT m =T c . If mitosis is measured in a population including quiescent cells, then the average proportion of cells in mitosis is lnð2ÞT m =T c multiplied by the growth fraction. Other asynchronous cell cycle distributions will change just the constant of proportionality in this expression, which is of consequence only if T m is of independent interest, as we will see below. A key problem is that T m (specifically defined as the period in which strong pHH3 marking occurs) is not known with any precision. Although believed to be on the order of 15–20 min (Cooper, 2000), the same proportional uncertainty must carry through to any estimate of T c using this formula directly. However it is widely assumed that T m is approximately constant for a given tissue type even if the overall cell cycle length varies, so the proportion of cells in mitosis is taken to be proportional (on average) to 1=T c , with a fixed constant of proportionality. Thus mitosis data may be used to measure variability in cell cycle length over time. If there is an undetermined level of quiescence or cell cycle length heterogeneity, we may use an alternative formulation in terms of proliferation rate. Consider the fact that the number of cells in mitosis is a direct measure of the increase in cell number via division over the succeeding time period of length T m (prior to consideration of entry and exit). Since T m is small, it follows that the proportion of cells in mitosis is a direct measure of the instantaneous proportional growth rate of the overall population (prior to entry and exit), regardless of growth fraction or other dynamics. However, whether working in terms of cell cycle length or proliferation rate, the unknown length of T m implies an unknown constant of proportionality. Determining absolute rather than relative values of cell cycle length or proliferation rate

72

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

requires an independent measure for at least one time point, such as the EdU analysis provides. A further issue in the use of pHH3 results is that the short length of T m translates to small proportionate numbers of pHH3 þ cells and relatively weak statistical power. In the best case the uncertainty can be modelled as a binomial process, however there may be a degree of local clonal synchronisation which reduces the effective sample size, as well as the possibility of systematic variability on the kidney or litter level. Use of pHH3 results therefore requires analysis of the reliability of any trend found. To experimentally measure the rate of death, a marker such as caspase may be used distinguish those cells that are undergoing apoptosis. Interpreting the proportion of apoptotic cells works similarly to mitotic cells; see Fig. 2 and discussion in Section 3.1. 2.4. Population growth and exit Since entry to the CM compartment has been excluded, the total cell population NðtÞ will satisfy dN ¼ ð∅ðtÞeðtÞÞN; dt where ∅ðtÞ is the proportional endogenous population growth rate (averaged over the total cell population including any quiescent cells) and eðtÞ is the proportional rate of exit. Thus R NðtÞ ¼ ke ∅ðtÞeðtÞdt ; for some constant of proportionality k. The function ∅ðtÞ is estimated using EdU and pHH3 results. Combined with CM population data over time, the rate of exit can be estimated over the course of development. Such an approach can be seen in (Bonhoeffer et al., 2000). Suppose that the population is known at times t 1 and t 2 ; say N 1 ¼ Nðt 1 Þ and N2 ¼ Nðt 2 Þ. Using the above formula for N, we have R t2 ∅ðtÞeðtÞdt: Nðt 2 Þ ¼ e t1 Nðt 1 Þ

Assuming that ∅ðtÞ and eðtÞ are constant over the interval ½t 1 ; t 2 , the exit rate is given by   Nðt 2 Þ =ðt 2 t 1 Þ: e ¼ ∅ ln Nðt 1 Þ Note that this depends only on the population ratio, hence for this calculation CM volume may be used as a proxy for population without any need to determine the constant of conversion between volume and cell count.

3. Results Here, the methods will be used to demonstrate two results in the cap mesenchyme, which yield insight into the turnover of cells in the compartment over the experimental time window. These results will show how the different experimental strategies may be combined within the mathematical framework described. 3.1. Data There are three pieces of experimental data on which the models will be based. For experimental and quantification methods, see Appendix 5.1. The first is the time course for the proportion of cap mesenchyme cells, as measured by being Six2 positive, that take up the cumulative S-phase marker over periods of up to of 16.5 h starting at 13.5 dpc. This is shown in Fig. 4A and will be used to estimate cell cycle and S-phase duration. The second is the total cap mesenchyme volume, as measured from 3D microscopy imaging of the Six2 positive cell volume, at daily intervals from 11.5 dpc to 17.5 dpc. This is shown in Fig. 4B and is used as a proxy for the number of cells within the cap mesenchyme at each time point. Thirdly, a static snapshot assessment of the proportion of cells in the cap mesenchyme that are undergoing mitosis at 24 h intervals from 11.5 dpc to 15.5 dpc is shown in Fig. 4C. A final consideration is apoptosis in the developing organ. As can be seen in Fig. 2D, there are very few cells at 13.5 dpc that

Fig. 4. Quantification results. In each experiment one litter was used for each time point, hence the standard deviation error bars are a measure of intra-litter variability only. Proportion of Six2 þ cells showing EdU incorporation, after course of EdU injections from 13.5 dpc. The horizontal scale is the hours since the initial injection (hpi). (B) Cap mesenchyme volume from OPT Six2 signal, measured in voxels. Each voxel is equal to a volume of 4.17 mm3.(C) Proportion of Six2 þ cells showing pHH3 incorporation. See Tables 4 and 5 (Appendix) for data values. (a) EdU uptake from 13.5 dpc, (b) Cap mesencyme volume per kidney and (c) Mitotic index in cap mesencyme.

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

are marked (red channel) as apoptotic. In particular, none are observed to be (green channel) cap mesenchyme cells. In multiple images containing a total of approximately 1000 cap mesenchyme cells, no caspase activity were observed within CM cells. This is consistent with other authors who have observed sparse apoptotic activity in mesenchymal or ureteric bud regions both pre and postnatally (Hartman et al., 2007). If it is assumed that apoptosis is detectable via the caspase-3 marker for an hour or more (Farfan et al., 2004), the proportion of cells undergoing apoptosis is negligible. Hence it will be assumed that there is no death during the experimental time window.

3.2. Calculation of cell cycle duration with sensitivity analysis In this section we analyse the results of the EdU experiment. Given the uncertainty regarding relative lengths of the cell cycle phases, timing of exit and possible heterogeneity of proliferation in the cap mesenchyme, we are faced with the task of fitting a model with a large number of possible parameters (some of which may have only a subtle effect on the EdU incorporation curve) to a relatively limited set of data. We take the approach of fitting several simplified models in order to provide a best estimate of the cell cycle length and proliferation rate, and also quantify the possible impact of undetermined parameters on these estimates.

73

As can be seen in Fig. 4A, the Six2 þ population exceeds 95% EdU incorporation at 10.5 hpi (10.5 h of exposure to EdU via intermittent injections, starting at 13.5 dpc), with a small increasing trend in the remainder of the time series. Direct determination of a saturation time is difficult to achieve with precision, with additional qualitative analysis of the image data at 12.5, 14.5 and 16.5 hpi (not shown) indicating no detectable long-term quiescence, but a small subpopulation of relatively slowly cycling cells which reach full EdU incorporation between 14.5 dpc and 16.5 dpc. Theoretical models as described in Section 2.2 were fitted to the data using the least-squares fitting function nls in the R language (R Core Team, 2013). Fitted curves and parameters are shown in Fig. 5 and Table 1, respectively. We first fitted single population models. A robust fit could not be obtained for the exponential model with T G2 þ m as a free parameter. With cell phase lengths constrained to be non-negative, the non-linear regression converged (giving T G2 þ m ¼ 0) but with standard errors greater than parameter values. The standard errors produced by the nonlinear least squares algorithm are an approximation based on a linearization of the parameter space around the fitted solution, and standard errors of this magnitude are not meaningful. To produce a robust fit, we repeated the regression with T G2 þ m ¼ 0 (Fig. 5A, solid line), giving T c ¼ 19:7 h. Although the length of the second gap phase and mitosis is known to be substantially shorter than the first gap phase, it is necessarily

Fig. 5. Selected models from Sections 2.2.1 and 2.2.3 fitted to EdU incorporation data of Fig. 4A. (A)Single population models, (Solid line: Exponential model with T G2 þ m ¼ 0, T C ¼ 19:7 h, Dashed line: Linear model, T C ¼ 17:2 h, Dotted line: Exponential model with T G2 þ m ¼ 2, T C ¼ 19:5 h). (B)Two population models, (Solid line: Exponential model with T G2 þ m ¼ 0, T C ¼ 19:8 h (average). Dashed line: Exponential model with T G2 þ m ¼ 0, subpopulations assumed to be of equal size, T C ¼ 15:4 h (average)).

Table 1 Parameters of fitted EdU models ( 7 standard error) shown in Fig. 5. Model

TC (h)

Ts (h)

TC1 (h)

TC2 (h)

Slow cycling cells (%)

Exponential model, T G2 þ m ¼ 0 h Linear model Exponential model, T G2 þ m ¼ 2 h Unconstrained two population model Two equal sized population model

19.7 7 2.0 17.2 7 1.6 19.5 7 2.1 19.8 7 1.7 15.4 7 3.1

8.0 7 1.5 5.9 7 1.1 8.0 7 1.7 8.0 7 0.8 5.7 7 1.8

18.5 7 1.7 12.2 7 3.2

25.9 7 6.1 20.8 7 2.3

237 32 50 (assumed)

74

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

greater than zero (Cooper, 2000). The length of T G2 þ m has been directly measured for various neural progenitor cells of comparable cell cycle length, using S-phase/mitosis double marking experiments. Lengths obtained were consistently 2 h or less (Takahashi et al., 1993, 1995; Cai et al., 1997; Calegari et al., 2005). Thus we also fitted the model with T G2 þ m ¼ 2 to test the sensitivity of the cell cycle length to this uncertainty, giving a noticeably worse fit (residual standard error increased from 0.037 to 0.043) with minimal change in T c or T s (Fig. 5A, dotted line). Although we lack independent evidence that T G2 þ m r 2 in the CM, higher values provide increasingly poor fits to the EdU data. As described in Section 2.2.1, we also test sensitivity to possible variation in the cell cycle distribution due to exit at mitosis, and to possible cell cycle synchrony of individual samples, by fitting the linear model. This gives a reduction in cell cycle length of only slightly over the margin of error (Fig. 5A, dashed line). The rapid expansion of the compartment suggests that the rate of exit from the compartment is significantly short of the 50% steady-state level consistent with the linear model, hence the uncertainty associated with the cell age distribution is minimal. The reduction in the fitted S-phase length is of similar magnitude but much more significant proportionally, so if the S-phase length is of interest then sensitivity to the cell cycle distribution assumptions is of much greater concern. The single population models provide a generally good fit to the data up to 10.5 hpi, but do not account for the small proportion of remaining EdU  cells at 12.5 hpi and 14.5 hpi, so we also considered two population models. As described in Section 2.2.3, we assume that S-phase length is common to the two cell populations in the CM, and T G2 þ m ¼ 0, to minimise the number of free parameters (TS, TC1, TC2 and a, the proportion of slower cycling cells). The resulting fit (Fig. 5B, solid line) resembles the corresponding single population model (Fig. 5A, solid line) with the addition of a small slow subpopulation. The fitted curves match almost exactly up until 10.5 hpi. Although the two population model gives a somewhat shorter cell cycle length for the main cell population compared to the corresponding single population model (also exponential with T G2 þ m ¼ 0), the average cell cycle length is almost identical (19.8 versus 19.7 h); recall that this average is a weighted harmonic mean of the subpopulation cell cycle lengths, giving the estimated doubling time of the entire cell population in the absence of entry and exit. This fit suggests that the presence of a small, slow cycling subpopulation can be ignored if we are interested in the overall proliferation rate, provided that the single population model provides a good fit to the major part of the time series. However, the two population fit, while providing a more convincing explanation of the data, was not entirely robust, with low confidence in the proliferation rate and particularly the size of the slower cycling subpopulation. To test the impact of this uncertainty about the relative sizes of the subpopulations, we repeated the regression with a ¼ 0:5 fixed (Fig. 5B, dashed line). The resulting fit is plausible, with full saturation reached between 14.5 dpc and 16.5 dpc as observed, and has only slightly greater residual standard error than the unconstrained two population fit (0.039 compared to 0.036). It represents a significantly different interpretation of the data, with the steeper slope of the early portion of the fitted curve corresponding to a higher proliferation rate. In summary, the one and two population models providing the best fit to the EdU data (exponential with T G2 þ m ¼ 0, solid lines in Fig. 5) give almost identical proliferation rates. There is little sensitivity to variation in length of T G2 þ m , or to timing of exit or synchrony of individual samples. The greatest uncertainty arises from cell cycle heterogeneity, as the available data suggests a small slower cycling subpopulation with little impact on overall proliferation, but cannot entirely exclude the possibility of two populations of similar size, associated with a significantly higher

aggregate proliferation rate. We proceed to the modelling of population dynamics using the single population exponential model with T G2 þ m ¼ 0, having an average cell cycle length of 19.7 h. 3.3. Predict change in proliferation rate over time We now use the cell cycle length estimate obtained from the EdU experiment to produce a model of population over time in the absence of quiescence or exit from the compartment (recall that entry has been excluded a priori). The average cell cycle length of 19.7 h ¼0.8211 days obtained from the EdU results corresponds to a proportional population growth rate of ln(2)/0.8211 ¼0.8441 days  1. Assuming a constant growth rate and extrapolating backwards and forwards in time from 13.5 dpc without exit, and continuing to use voxels as a proxy for population size, we get a population size at t days dpc of NðtÞ ¼ 4488798e0:8441ðt13:5Þ ¼ 50:5e0:8441t This model will be referred to as the constant model. The constant model assumes that the proliferation rate within the CM is invariant. Using the mitosis data (Fig. 4C) we investigate the possibility that it may instead systematically vary with time. The mitotic index measurements are highly variable, but there is evidence for a decreasing trend. We obtain a linear regression coefficient of 1:62  103 7 4:6  104 , with a p-value of 0.04 (Fig. 6), indicating a slowdown in proliferation over time but with high uncertainty. Note that as stated in Section 2.2, the mitotic index is a direct proxy for overall proliferation, and this trend can be used regardless of whether the observed slowdown is due to increasing cell cycle length or to cells exiting the cell cycle and entering a quiescent state (while remaining within the CM) after 13.5 dpc. However, the absence of quiescent cells at 13.5 dpc implies that quiescence cannot explain the trend of decreasing mitosis over the entire timecourse from 11.5 dpc to 15.5 dpc, so increased cell cycle length is the simplest explanation for this trend. Extrapolating from the proliferation result obtained from the EdU experiment at 13.5 dpc (single population model with T G2 þ m ¼ 0), we obtain a proliferation rate over time of ∅ðtÞ ¼ 2:30470:1082t, giving the variable model R 2 NðtÞ ¼ ke ∅ðtÞdt ¼ 0:00264e2:3047t0:0541t :

Fig. 6. Linear trend line fitted to average mitotic index.

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

3.4. Model comparison The constant and variable models can be used to find evidence for proliferation change and differentiation, by comparing the models of population size predicted to the observed values. The difference between the model and the observed values provides information about any unaccounted-for behaviours in the compartment that would affect the population size over time. In Section 3.3 the constant and variable models were obtained by fitting to the experimental proliferation data in the cap mesenchyme, and the observed CM volume at 13.5 dpc. Note that these models do not account for any potential exit. In Fig. 7 the models are compared to the observed population values.

3.4.1. Evaluation of possible increase of cell cycle length between 11.5 dpc and 14 dpc The constant model provides a good fit to the data at 12 dpc, 12.5 dpc, 13.5 dpc and 14 dpc, but predicts a larger population at 11.5 dpc than the observed values (Fig. 7B). The variable model provides a better fit to the data at 11.5 dpc, while giving a similarly good fit at the other time points up to 14 dpc. Thus the 11.5 dpc population data supports the variable model, with a shorter cell cycle before 13.5 dpc. However, it should be noted that the cap mesenchyme at 11.5 dpc is much less compact than at later time points, which may indicate that this data is less reliable. Thus the population data up to 14 dpc is consistent with a trend of increasing cell cycle length, but is not definitive.

3.4.2. Evidence for differentiation or quiescence from 15.5 dpc to 17.5 dpc The constant model predicts a much larger population from 15.5 dpc to 17.5 dpc than the observed values (Fig. 7A), 3.6 times as large by 17.5 dpc. This might be explained by either reduced

75

proliferation in the cap, due to increased cell cycle length, the establishment of a substantive quiescent fraction within the cap mesenchyme, or by exit from the Six2 positive compartment via differentiation. As noted in Section 3.1 above, there is little evidence of significant apoptosis within the developing cap mesenchyme, and so any exit is assumed to be by differentiation. The variable model based on the mitosis data does predict an increasing length of cell cycle as described in 3.3, and the data prior to 13.5 dpc is consistent with this model, as noted above. However, while the model does indeed predict cap volumes closer to the observed values than the constant model for 15.5–17.5 dpc (Fig. 7A), there is a significant and consistently increasing difference between the variable model and the observed values, with the predicted population 50% larger than that observed by 17.5 dpc. So even if we assume a significant ongoing decline in proliferation rate, we see a discrepancy much larger than can be explained by population measurement error. Hence, variation in cell cycle length is not sufficient to explain the divergence from exponential growth. The two remaining possibilities are a significant quiescent population or significant differentiation. We examine and quantify each of these cases. First, suppose that there is no differentiation, and that the difference between the models and the experimental data is entirely due to a proportion of the cells at each time point becoming quiescent (assuming exit from the cell cycle but not from the compartment). For each day between 13.5 dpc and 17.5 dpc, we calculate the proportion of cells that must be quiescent in order to account for the difference between the predictions of the variable or constant model on the one hand and the experimentally determined cap mesenchyme volumes on the other. The proliferating population is assumed to follow either the constant or variable model. We give the calculation for the constant model below; the variable model follows similarly. For a given day, define the population volume to be P0 ¼ QþA0, where Q is the quiescent population at that time point.

Fig. 7. Comparison of experimentally quantified (see Fig. 4B) and predicted compartment volumes (voxels) across developmental time extrapolated based upon the observed volume and calculated cell cycle length Tc at 13.5 dpc. Predicted compartment volumes are calculated using the constant model (proliferation unchanged across time) or the variable model (proliferation varies according to observed changes in mitotic index across time), as given in Section 3.3. (A) Across full time series 11.5–17.5 dpc. (B) Detail, 11.5–14 dpc only.

76

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

The subsequent day, in the constant model, should then have population P1 ¼Qþ A0e0.8411. Given the experimentally measured values for P0 and P1 (Fig. 4 and Table 4 below), it is straightforward to solve for Q and A0, and then determine Q/P0, the fraction of cells that are quiescent on the given day. The result of these calculations using the constant and variable models on days 13.5 dpc to 16.5 dpc are given in Table 2. It can be seen that the predicted proportion of quiescent cells, under the assumption of no exit from the compartment, are generally quite high for both models with the minimum quiescent proportion for either model of 11%. However, the S-phase EdU marker experiments given in Fig. 5 for 13.5 dpc show that almost all cells are passing through S-phase and hence that the quiescent population is of negligible size at this time point. This suggests that the primary mode for the reduction in growth rate in the cap mesenchyme at later time points is through differentiation and exit from the compartment rather than quiescence, although a combination of factors cannot be ruled out. Suppose instead, that there is no quiescence in the CM, but that the difference between the models and the experimental data is in fact due to exit by differentiation, a process known to be occurring by 12.5 dpc (Park et al., 2007). The rate of differentiation (or of total exit, including cell death also) is calculated as described in Section 2.4. For each day between 13.5 dpc and 17.5 dpc, we assume a fixed exit rate x for that day. We fit the exponential function NðtÞ ¼ Aeð0:8441xÞt (for the constant model), to the data points for that day and its successor, to solve for A and the exit rate x for that day. For instance, fitting N(t) to the experimentally measured volumes for 14.5 dpc and 15.5 dpc, (uniquely) determines an exit rate of x ¼0.218. Similarly, exit rates for the variable model using NðtÞ ¼ Aeð2:30470:0541txÞt can be uniquely determined. In Table 3, the list of exit rates, and hence the predicted proportion of cells differentiating at each time point, are given. For the constant model, up to 50% of cells in a given day are predicted to differentiate. The proportion of cells predicted to differentiate is lower for the variable model, with an exit rate of 0.231 between 15.5 dpc and 16.5 dpc. Using these fitted exit rates and the population sizes, we can estimate the daily volume of exit from the CM. Fig. 8 shows population and exit from 13.5 dpc for the variable model. The constant model, without the assumption of increasing cell cycle length, produces much higher estimates of exit. Reliable quantification of the amount requires more precise estimation of proliferation later in the course of development, to determine the correct model. But in either case, we have shown there is strong evidence for significant exit from the compartment.

Table 2 Predicted quiescence in the cap mesenchyme compartment as percentage of population per day. dpc Percentage of in quiescence Percentage of in quiescence

cells Constant model cells Variable model

13.5–14.5 20.3%

14.5–15.5 34.4%

15.5–16.5 69.2%

16.5–17.5 63.2%

12.1%

11.0%

47.2%

17.7%

Table 3 Predicted exit rate from the cap mesenchyme compartment as fraction of population per day. The exit rate corresponds to the proportion of cells in that day that are exiting. dpc Constant model exit rate Variable model exit rate

13.5–14.5 0.123 0.069

14.5–15.5 0.218 0.056

15.5–16.5 0.501 0.231

16.5–17.5 0.447 0.068

Fig. 8. Observed CM populations from 13.5 dpc, in millions of voxels, and the modelled exit over the previous day on the same scale. Assumes variable model and no quiescence.

4. Conclusion A novel method for quantifying and modelling cell cycle length, population heterogeneity and developmental exit via differentiation within an organ compartment has been demonstrated. A key element has been to use a model of the cumulative S-phase marker uptake over time that is appropriate to the assumptions and timescales typically occurring in the rapidly changing environment of a developing organ and to be able to estimate the affect of potential confounding factors such as cell cycle synchrony. The cell cycle and S-phase durations may then be estimated from the model. These then form the foundation for an exponential growth model of cell proliferation in the organ. Experimentally measured deviations in cell population size from this model then need to be explained by mechanisms such as cell differentiation, apoptosis, variation in cell cycle length, quiescence and compartment entry. Here we demonstrated how this might be done using mitosis and apoptosis data. This led to two key results; that the cell cycle increases in length from 11.5 dpc to 13.5 dpc and that there is good evidence for significant differentiation from 15.5 dpc to 17.5 dpc. It should be noted that the aim of the paper has not been to conclusively prove these results. Rather, the aim of the paper has been to show how each of these mechanisms for cell population change can be measured, modelled and accounted for in order to understand the dynamic environment of a developing organ, and to outline the requirements and priorities for a more definitive analysis, as follows. This method can also afford a better resolution of processes contributing to an abnormal state that the more common approach of relative comparisons of simpler measurements such as proportion of mitosing cells, the latter demonstrated here to represent a very shallow resolution tool. First, due to the variability of mitosis measurements, cumulative S-phase marker uptake experiments could be performed on multiple days to more reliably quantify cell cycle length and proliferation, or quiescence, over the course of development.

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

This could include an experiment performed at the latest practical stage, in order to quantify or exclude quiescence in the compartment as reliably as possible. However, these are costly experiments in time and materials and so it is highly useful to have the strong predictive mathematical evidence, such as the models we provide here, before performing them. Regarding the cumulative S-phase marker uptake experiment itself, special care should be taken if substantial heterogeneity of proliferation at a given time is considered biologically plausible. For example, if the compartment comprises two cell subpopulations of comparable size and significantly different cell cycle lengths. In this heterogeneous case substantially greater precision is required in the time series to provide a robust fit. While a single population model can be reasonably applied given only the initial incorporation level and the saturation time, a two population model requires accurate measurements throughout the time course, with the estimate of overall proliferation rate heavily dependent on the slope prior to the saturation of the faster cycling subpopulation, as demonstrated in Section 3.2. The minimal requirement for quantifying two subpopulations is two measurements prior to the saturation of the faster subpopulation, and a further two measurements between the saturation of the faster population and the saturation of the slower population. For a compartment with similar characteristics to the CM, we recommend using time points at a maximum spacing of 4 h from initial injection until saturation has been established, with multiple litters at each time point if possible. Variation in the length of the second gap phase has a fairly subtle effect on the predicted incorporation curve. As a result it is difficult to determine by fitting to the observed incorporation, but conversely it has little impact on the cell cycle length estimate. If desired, this parameter may be estimated more directly using a cumulative S-phase and mitosis duel marker experiment, see for example (Calegari et al., 2005). We also demonstrated that the uncertainty due to timing of exit or due to cell cycle synchrony in individual samples is a minor concern in the CM (or a compartment with similar characteristics), while giving a method to calculate a bound on this uncertainty. While our approach has been demonstrated on the cap mesenchyme compartment within the developing mouse kidney, and the above recommendations developed in that context, it should be applicable to a wide range of progenitor populations or cellular states in developmental biology.

Author contributions JL and DJM developed the mathematical models of the paper and substantially wrote the manuscript. NAH, MHL and ANC contributed to the intellectual design of the research and also wrote substantial parts and contributed to the figures. ANC designed and performed the bulk of the experimental work together with ALJ. All authors read and approved the final manuscript before submission.

Acknowledgements ML is an NHMRC Senior Principal Research Fellow. This work was supported by the National Health and Medical Research Council (APP1002748) and the Human Frontiers Science Program (RGP0039/2011). Confocal microscopy and optical projection tomography were performed at the Australian Cancer Research Foundation Cancer Biology Imaging Facility. We wish to acknowledge the facility manager, James Springfield, for designing and constructing the OPT hardware and software.

77

Appendix Experimental and quantification methods Mouse strains and immunofluorescence Six2-TGC BAC mice (Kobayashi et al., 2008)were mated with CD1 mice, with noon of the day on which the mating plug was observed designated 0.5 dpc. Embryos that inherit the Six2-TGC BAC transgene from these matings express green fluorescent protein (GFP) at the sites of native Six2 expression (Kobayashi et al., 2008). Therefore Six2þ cells can be identified with an antibody to GFP protein. Whole mount immunofluorescence for confocal was performed on isolated kidneys as described (Combes et al., 2009) but samples were dehydrated in a MeOH series and cleared in a 1:2benzyl alcohol: benzyl benzoate solution prior to imaging in a glass bottom dish (Mattek P35G-1.5-14-C). Sample preparation for OPT was as described (Rumballe et al., 2011). The following primary antibodies were used at a dilution of 1:300: anti-GFP (Abcam Ab13970), antiCalbindin (Sigma c-9848), anti-Phh3 (Millipore 06–570), antiCaspase3 (CST 9661S). Secondary Alexa Fluor-conjugated antibodies and DAPI (Invitrogen) were used at 1:500. Optical projection tomography to estimate Six2 volume Six2 marks the cap mesenchyme population although it is also expressed at lower levels in early nephrons (Rumballe et al., 2011). For optical projection tomography (OPT), multiple samples were collected from between one and three litters and data was averaged. Number of samples per age is as follows: 11.5, 7; 12.0, 2; 12.5, 10; 13.5, 11; 14, 5; 14.5, 6; 15.5, 9; 16.5, 3; 17.5, 3. For OPT, the 11.5 dpc to 15.5 dpc samples were imaged with a  1 objective and  4 zoom; the16.5 dpc samples were imaged with a  1 objective,  3 zoom; and the 17.5 dpc samples with a  1 objective,  2 zoom. In the Six2 OPT imaging, there is a low cellular resolution but a full kidney in view. Because the nuclei are packed in tightly, it is assumed that the nuclei, as a combined mass, denote regions that are entirely Six2þ . To determine the amount of Six2 present, total volume of Six2 in the whole kidney is quantified. From an OPT image stack, ImageJ (Abramoff et al., 2004) was used to segment the Six2þ region. A median filter of radius 1 was performed, then for the samples taken with a  1 objective and  4 zoom, background subtraction of 20 and a threshold 35 was used to define the positive regions on each slice. For the  1 objective and  3 zoom a background subtraction of 12 and threshold of 35 was used and for the  1 objective,  2 zoom samples, background subtraction of 8 and a threshold of 35 was used. The number of pixels in each slice was then counted. Since the X, Y and Z scales are the same, each pixel is a 1  1  1 voxel and then total voxel count signifies the volume of cap mesenchyme in that particular kidney. To use the volume as a proxy for cell count, it is required that the average cell size and density does not change significantly across time. In the Six2 cap mesenchyme data, this is assumed. The voxel counts are also corrected with respect to the change in objectives as the scale in the imaging changes. This is done by multiplying through by the difference between the scale for the  1 objective,  2 zoom, compared with the  3 and  2 zooms. In this way, all the data is in terms of the units of scale of the  1 objective,  4 zoom. For instance, a 1  1  1 voxel in the  3 zoom imaging is the same size as a (4/3)3voxel in the  4 zoom imaging. Such a pixel will have an volume of 64/27 in the  4 zoom scale. At the  4 zoom, the side length of a pixel is 1.61 mm, giving a volume per voxel of 4.17 mm3. EdU uptake in Six2 population An amount of 1 mg EdU was delivered via intraperitoneal injection at 9am when mice were 13.5 dpc. The initial injection

78

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

Table 4 Average population volume and mitotic index over time for the cap mesenchyme. See Fig. 4B and C for plots. dpc # kidneys analysed for Six2 volume Mean Six2 volume (voxels) Standard deviation # kidneys analysed for mitosis Mean mitotic index Standard deviation

11.5 7 420862 265539 2 0.0156 0.0090

12 2 1170317 127990 – – –

12.5 10 1833695 785477 3 0.0151 0.0023

time was designated 0 h and 1 mg boost injections were delivered every 4 h. The 4 h time point was selected based on experiments demonstrating a loss of EdU staining in the cap mesenchyme at greater than 5 h (data not shown).A series of litters were taken at times (litter size) 0.5 h (5), 4.5 h (4), 8.5 h (5), 10.5 h (3), 12.5 h (5), 14.5 h (3), 16.5 h (5). For the litters taken at 10.5 h and 14.5 h, an additional injection was delivered at 10 h and 14 h respectively, ensuring that there is always an injection 30 min prior to harvesting. This was based on experiments consistently showing a plateau or small decline in staining at these time points when no additional injection was delivered, indicating that the EdU persists in the system for significantly less than 2 h (data not shown). We used 30 min after injection as our first time point to ensure that EdU had sufficient time to reach the embryos. EdU, and GFP (Six2GFP/cap mesenchyme) signal was imaged on an inverted Zeiss LSM 510 Meta confocal microscope. Samples were imaged in nonoverlapping 2.03 μm optical slices for up to 120 μm. The image analysis program Imaris was used to measure EdU levels within Six2 positive cells. Briefly, workflow and analysis was as follows: EdU data was manually analysed to determine an intensity threshold cutoff that clearly separated EdU signal from background, the EdU signal was then binarised. Six2-GFP signal was used to surfacerender Six2 expressing nuclei, which were subsequently filtered by Six2 intensity (to exclude nephron-based Six2-GFP expression) and by cell size (250–1500 um3, to exclude clusters of cells). The binarised EdU was then measured within each remaining distinct Six2-GFP nucleus, and the nucleus was classified as EdU þ if greater than 20% of the pixels within that volume was occupied by EdU signal. This process was repeated with multiple samples for a time point and data expressed as an average % of the population where Six2-GFP nuclei were considered EdU þ ve.

Mitosis in the Six2 population Phospho histone H3 (pHH3) is a marker that is enriched in condensed chromosomes during the mitosis (M) phase of the cell cycle. Immunofluorescence for Phh3 and Six2-GFP was performed on isolated kidneys to determine the proportion of cap mesenchyme cells undergoing mitosis at one time. Litters were taken at times (litter sizes) 11.5 dpc (2), 12.5 dpc (3), 13.5 dpc (4), 14.5 dpc (3) and 15.5 dpc (3). Confocal imaging was as described above. Imaris was used to count the number of Six2-GFP nuclei in a field of view, and to count which of those cells were also expressing Phh3. To identify Six2-GFP nuclei a manual threshold was set to identify Six2-GFP signal above background, and nuclei were identified using a spot counting algorithm using a size guide of 3 μm. A 3D colocalisation was performed to identify Six2-GFP/ Phh3 double-positive cells which were then counted. Data was expressed as an average % of double positive nuclei compared to the total.

Data tables See Tables 4 and 5.

13.5 11 4488798 1188328 4 0.0111 0.0056

14 5 6602763 596699 – – –

14.5 6 9235365 1757937 3 0.0127 0.0019

15.5 9 17267749 3830628 3 0.0087 0.0014

16.5 3 24330325 1480721 – – –

17.5 3 36193224 6521543 – – –

Table 5 Proportion of Six2 þ cells showing EdU incorporation, after course of EdU injections from 13.5 dpc. See Fig. 4A for plotted data. hpi # kidneys analysed Total # Six2 þ cells EdU proportion Standard deviation

0.5 5 12784 0.343 0.018

4.5 4 13128 0.660 0.047

8.5 5 10741 0.800 0.067

10.5 3 11709 0.955 0.013

12.5 5 12592 0.955 0.025

14.5 3 9927 0.971 0.028

16.5 5 4328 0.991 0.005

References Abramoff, M.D., Magelhaes, P.J., Ram, S.J., 2004. Image processing with ImageJ. Biophotonics International 11 (7), 36–42. Alexiades, M.R., Cepko, C., 1996. Quantitative analysis of proliferation and cell cycle length during development of the rat retina. Developmental Dynamics 205 (3), 293–307. Baker, C.T., et al., 1998. Modelling and analysis of time-lags in some basic patterns of cell proliferation. Journal of Mathematical Biology 37 (4), 341–371. Bertuzzi, A., et al., 1997. Steel's potential doubling time and its estimation in cell populations affected by nonuniform cell loss. Mathematical Biosciences 143 (2), 61–89. Bonhoeffer, S., et al., 2000. Quantification of cell turnover kinetics using 5-bromo2′-deoxyuridine. The Journal of Immunology 164 (10), 5049–5054. Boyle, S., et al., 2008. Fate mapping using Cited1-CreERT2 mice demonstrates that the cap mesenchyme contains self-renewing progenitor cells and gives rise exclusively to nephronic epithelia. Developmental Biology 313 (1), 234–245. Cai, L., Hayes, N.L., Nowakowski, R.S., 1997. Local homogeneity of cell cycle length in developing mouse cortex. The Journal of Neuroscience 17 (6), 2079–2087. Calegari, F., et al., 2005. Selective lengthening of the cell cycle in the neurogenic subpopulation of neural progenitor cells during mouse brain development. The Journal of Neuroscience 25 (28), 6533–6538. Caviness Jr., V.S., Takahashi, T., 1995. Proliferative events in the cerebral ventricular zone. Brain & Development 17 (3), 159–163. Charvet, C.J., Striedter, G.F., 2008. Developmental species differences in brain cell cycle rates between northern bobwhite quail (Colinus virginianus) and parakeets (Melopsittacus undulatus): implications for mosaic brain evolution. Brain, Behavior and Evolution 72 (4), 295–306. Combes, A.N., et al., 2009. Three-dimensional visualization of testis cord morphogenesis, a novel tubulogenic mechanism in development. Developmental Dynamics 238 (5), 1033–1041. Cooper, G.M., 2000. The Cell, second ed. Sunderland (MA): Sinauer Associates. Farfan, A., et al., 2004. Multiplexing homogeneous cell-based assays. Cell Notes 10, 15–17. Georgas, K., et al., 2009. Analysis of early nephron patterning reveals a role for distal RV proliferation in fusion to the ureteric tip via a cap mesenchymederived connecting segment. Developmental Biology 332 (2), 273–286. Hartman, H.A., Lai, H.L., Patterson, L.T., 2007. Cessation of renal morphogenesis in mice. Developmental Biology 310 (2), 379–387. Hendry, C., et al., 2011. Defining and redefining the nephron progenitor population. Pediatric Nephrology 26 (9), 1395–1406. Kobayashi, A., et al., 2008. Six2 defines and regulates a multipotent self-renewing nephron progenitor population throughout mammalian kidney development. Cell Stem Cell 3 (2), 169–181. Kornack, D.R., Rakic, P., 1998. Changes in cell-cycle kinetics during the development and evolution of primate neocortex. Proceedings of the National Academy of Sciences of the United States of America 95 (3), 1242–1246. Mugford, J.W., et al., 2009. High-resolution gene expression analysis of the developing mouse kidney defines novel cellular compartments within the nephron progenitor population. Developmental Biology 333 (2), 312–323. Nowakowski, R.S., Lewin, S.B., Miller, M.W., 1989. Bromodeoxyuridine immunohistochemical determination of the lengths of the cell cycle and the DNA-synthetic phase for an anatomically defined population. Journal of Neurocytology 18 (3), 311–318. Park, J.S., Valerius, M.T., McMahon, A.P., 2007. Wnt/beta-catenin signaling regulates nephron induction during mouse kidney development. Development 134 (13), 2533–2539. R Core Team, R: A Language and Environment for Statistical Computing. 2013, R Foundation for Statistical Computing: Vienna, Austria.

J. Lefevre et al. / Journal of Theoretical Biology 338 (2013) 66–79

Rumballe, B.A., et al., 2011. Nephron formation adopts a novel spatial topology at cessation of nephrogenesis. Developmental Biology. 360 (1), 110–122. Steel, G.G., 1967. Cell loss as a factor in the growth rate of human tumours. European Journal of Cancer (1), 381–387. Steel, G.G., 1977. Growth Kinetics of Tumors: Cell Population Kinetics in Relation to the Growth and Treatment of Cancer. Clarendon Press, Oxford, UK. Steel, G.G., Bensted, J.P.M., 1965. In vitro studies of cell proliferation in tumours I: Critical appraisal of methods and theoretical considerations. European Journal of Cancer (1), 275–279.

79

Takahashi, T., Nowakowski, R.S., Caviness Jr., V.S., 1993. Cell cycle parameters and patterns of nuclear movement in the neocortical proliferative zone of the fetal mouse. The Journal of Neuroscience 13 (2), 820–833. Takahashi, T., Nowakowski, R.S., Caviness Jr., V.S., 1995. The cell cycle of the pseudostratified ventricular epithelium of the embryonic murine cerebral wall. The Journal of Neuroscience 15 (9), 6046–6057.