An Unbiased Analysis of Candidate Mechanisms for the Regulation of Drosophila Wing Disc Growth

The control of organ size presents a fundamental open problem in biology. A declining growth rate is observed in all studied higher animals, and the growth limiting mechanism may therefore be evolutionary conserved. Most studies of organ growth control have been carried out in Drosophila imaginal discs. We have previously shown that the area growth rate in the Drosophila eye primordium declines inversely proportional to the increase in its area, which is consistent with a dilution mechanism for growth control. Here, we show that a dilution mechanism cannot explain growth control in the Drosophila wing disc. We computationally evaluate a range of alternative candidate mechanisms and show that the experimental data can be best explained by a biphasic growth law. However, also logistic growth and an exponentially declining growth rate fit the data very well. The three growth laws correspond to fundamentally different growth mechanisms that we discuss. Since, as we show, a fit to the available experimental growth kinetics is insufficient to define the underlying mechanism of growth control, future experimental studies must focus on the molecular mechanisms to define the mechanism of growth control.

Here, we use published quantitative growth data to carry out an unbiased evaluation of alternative candidate growth mechanisms. Based on the quantitative analysis, we rule out growth control by dilution in the wing disc. Evaluating alternative candidate growth mechanisms, we find that biphasic exponential growth best fits the wing disc data. An exponential decline of the growth rate with developmental time and logistic growth are also consistent with the data. These three growth laws correspond to fundamentally different growth mechanisms that we discuss. Given the wide range of growth laws that are consistent with the data, we emphasize the need to explore a wider range of growth limiting mechanisms and to carefully check the consistency of any proposed mechanism with additional experimental data -the ability to recapitulate the growth kinetics alone provides insufficient support for any mechanism.

Results and Discussion
Dilution of a cytokine cannot explain the Wing Disc Growth Kinetics. The growth of the apical area, A, of Drosophila imaginal discs over developmental time, t, can be described mathematically by Here, k refers to the area growth rate. We have previously shown that the declining area growth rate in the eye disc can be described by an area-dependent decline 34 of the form Here A 0 is the initial area and k 0 the maximal growth rate. We sought to test whether this mechanism would also apply to the wing disc. To this end, we used two independent, published datasets for wing disc area growth 17 1A) to determine the area growth rate, k (Fig. 1B), and to fit the model (Equations 1 and 2) to these datasets (Fig. 1C,D). The poor fit to the data rejects this mechanism for growth control in the wing disc (Fig. 1C,D, Figs S1A-D and S2A-D).
Unbiased Analysis of Candidate Growth Laws for Wing Disc Growth. We next tested which other candidate growth laws would allow us to reproduce the wing disc data (Fig. 2). To this end, we simulated the growth model in Eq. (1) with the growth laws that we had previously considered as candidates for growth control in the Drosophila eye disc 34 . In addition, we considered biphasic growth, which had previously been shown to describe wing disc growth 2 , and logistic growth, which had previously been shown to describe the regeneration kinetics of newt limbs 14 . In summary, in the simplest model we assumed a constant growth rate (CST, Fig. 2A In biphasic growth (BPH), the value of k 0 changes once, thus giving rise to two phases of different exponential growth (biphasic, BPH, Fig. 2C,D). The switch point from phase I to II was defined as the one providing the lowest deviation from the data (Fig. S3). In addition, we considered models with a continuously declining growth rate, either as exponential decaying growth rate (EXP, Fig. 2E,F) of the form Finally, the growth rate could depend on the total wing disc area. Besides dilution (Equation 2, Fig. 1C,D), logistic growth (LOG, Fig. 2I,J) can result in an area-dependent growth rate We note that the models with a constant or an area-dependent growth rate have two parameters, the biphasic exponential model (BHP) has four parameters, while all other models (EXP, POW, LOG) have three independent parameters. A full list for the parameter values can be found in the supplementary information (Tables S1 & S2).
To fit the models to the data, we defined the residuals in two independent ways: either by normalizing with respect to the standard error (dashed lines in Figs Table S3) and coefficient of determination (R 2 , Table S4) were used. We note that fitting the logarithmic data provided the smaller deviation and better coefficient of determination (R 2 ) in all cases (Tables S3 and S4). Using these measures, the biphasic model (BHP) performs best for almost all weightings and data sets (Fig. 2K). Interestingly, also the reported change in cell density 17 correlates with the predicted switch point. Thus, the cell density increases until the inferred switch point (vertical lines) and subsequently stabilizes (Fig. 2L). Logarithmic growth (red) and to a slightly lesser extent an exponentially declining growth rate (yellow) also fit the data very well (Fig. 2K), and the difference to the biphasic model is only minor (Tables S3 & S4). The power law model fits the data slightly, but consistently, worse as judged by the R 2 value (Table S4) and the RSS (Table S3).
The models with a constant (CST) or an area-dependent growth rate both have only two free parameters, and thus at least one parameter less than the models discussed above (BHP, EXP, POW, LOG). The model with the area-dependent growth rate performs considerably worse than the model with a constant (CST) growth rate, and can thus be rejected (Fig. 2K, Tables S3 & S4). To compare the model with a constant (CST) growth rate to the other models, the different number of free parameters needs to be taken into account. Unlike other commonly used model selection criterions, such as the Bayesian information criterion (BIC) or the Akaike information criterion (AIC), the F-test provides a p-value on the null hypothesis that the better fit can be solely explained by the increased number of parameters. Its use is, however, restricted to the case where the simpler model is nested within the more complex model, limiting it to the comparison of the biphasic (BHP) and single exponential model (CST). Conducting an F-test, we find that for both data sets and weightings, the biphasic exponential model fits the data significantly better than the model with a constant growth rate (Table S5).

Candidate Mechanisms for Growth Control in the Drosophila Wing Disc.
While the biphasic growth law, logistic growth, and an exponentially declining growth rate all fit the data well, they point to very different underlying mechanisms. A biphasic growth law would require a sudden change in cell growth, potentially because of cell differentiation. It has been argued that larvae can monitor their size and trigger the switch to lower growth rates when reaching a 'critical size' 2 . We note that in contrast to models with continuously declining growth rates, the biphasic growth law would still require an additional mechanism to ultimately stop growth. Otherwise, growth would continue with the speed from the second phase. Growth in the models with a continuously declining growth rate will eventually be so close to zero that the expansion becomes negligible.
An exponential decline in the growth rate could be achieved by the Dpp-dependent mechanism proposed by Gonzalez-Gaitan and co-workers 17,28 , or if a growth factor to which the system responded linearly, was degraded at a constant rate δ , i.e.
However, in both cases there are fundamental problems such as the reliable read-out at low concentrations and the robustness of the mechanism to changes in the total developmental time as observed in grafting experiments 10,11 . Thus, as the exponentially declining growth rate, k (Equation 4), declines 3-fold and 5-fold in the two different datasets from early to late stages 17,35 , cells would need to be able to linearly respond to 3-5-fold changes in this growth-controlling factor. Recent experiments indeed question a central role of the Dpp gradient in determining final disc size [29][30][31] .
Logistic growth requires a mechanism to set a final size and to reduce the growth rate as this final size is approached. The intercalation model results in logistic growth 14 . Here, the growth rate is postulated to be proportional to the positional difference between neighbouring cells, which is reduced as the structure grows out. How such positional differences would be measured by cells is not known. While loss of the cell-adhesion mediating cadherin Fat has been shown to enhance growth in the medial part of the wing disc, it has remained unclear whether the Fat-dependent growth limitation is central to growth termination 33,[36][37][38] . Beyond its molecular implementation, an intercalation mechanism has a number of conceptual limitations. Thus, while the intercalation mechanism can explain restoration of missing positional values during regeneration and a logistic growth law would be robust to changes in developmental speed, it remains unclear how the positional identity would scale when tissues grow to different final sizes, for instance because of differences in nutrients available 2 . Finally, while it has been argued that larvae can sense a 'critical size' 2 , it remains unclear how wing discs would achieve this on the molecular level in a disc-autonomous way and how disc size could then vary in response to changes in external conditions (nutrients, temperature etc.).
Scientific RepoRts | 6:39228 | DOI: 10.1038/srep39228 In conclusion, we can rule out an area-dependent dilution mechanism for growth control in the wing disc. Rather the responsible mechanism must give rise to biphasic growth, logistic growth, or an exponentially declining growth rate. Given that several growth laws match the measured growth curves, it needs to be stressed that the reproduction of the growth kinetics alone is insufficient evidence for any proposed mechanism and mechanistic proof needs to be provided. At the same time, there could be other growth laws that also reproduce the data and that we have not yet identified. It will be interesting whether the same mechanism limits growth in all organs and species or whether distinct mechanisms have evolved 17,28,34,39,40 .

Materials and Methods
Software. All simulations and analysis were done using Matlab R2016A (The MathWorks, Natick, MA, USA) or the free software environment R in version 3.2.4. 41 .
Estimation of area growth rate. The area growth rate 8) in Eq. (1) can be determined by estimating the slope dA dt at all time points and diving this slope by the measured area A at each timepoint. To estimate the slope dA dt from the data without making any assumption about the (true) function underlying the growth dynamics, we used cubic smoothing splines as provided by the smooth.spline function in R 41 . As the exact shape of the decline of the area growth rate k depends on the spline fits, we used two different fits for each data set with a varying smoothing parameter. The resulting spline fits are shown in Fig. 1A, and the estimated area growth rate k(t) in Fig. 1B. Model Fitting. The two available independent studies do not provide the variance for all data points 17,35 . We therefore quantified the deviation of the model from the data in form of the residual R at time point k as 10 10 or by weighting the difference by the given standard error (SE) as where D k is the mean, SE k its standard error, and f(x,t k ,p) refers to the model value at time point k. We took the logarithm of the data for the first measure since the data is spread over several orders of magnitude and larger values would therefore be weighted much higher otherwise. In the data set from Wartlick et al. no SE was given for the time point at approximately 80 hours and for the last time point (SE = 0) 17 . We therefore excluded these two time points from the analysis with Eq. 10. The model parameters were then determined by solving the minimization where N is the number of distinct time points. This was done independently for both definitions of the residuals using the trust-region reflective algorithm as implemented in lsqnonlin (Matlab R2016A, MathWorks, Natick, MA, USA), resulting in two parameter sets for each model. where RSS j is the residual sum of squares of model j calculated as explained above (Equations 9 or 10) and SST is the total sum of squares.
To evaluate whether there is a statistical significance in the difference how the models fit the data, the F-test was used. The F-test can only be used to compare a simple model which is nested within a more complex one. Thus, here it can only be used to compare the biphasic exponential model to the single exponential model, as the latter one is nested within the first one. The F-statistic was calculated as