FGF2 modulates simultaneously the mode, the rate of division and the growth fraction in cultures of radial glia

ABSTRACT Radial glial progenitors in the mammalian developing neocortex have been shown to follow a deterministic differentiation program restricted to an asymmetric-only mode of division. This feature seems incompatible with their well-known ability to increase in number when cultured in vitro, driven by fibroblast growth factor 2 and other mitogenic signals. The changes in their differentiation dynamics that allow this transition from in vivo asymmetric-only division mode to an in vitro self-renewing culture have not been fully characterized. Here, we combine experiments of radial glia cultures with numerical models and a branching process theoretical formalism to show that fibroblast growth factor 2 has a triple effect by simultaneously increasing the growth fraction, promoting symmetric divisions and shortening the length of the cell cycle. These combined effects partner to establish and sustain a pool of rapidly proliferating radial glial progenitors in vitro. We also show that, in conditions of variable proliferation dynamics, the branching process tool outperforms other commonly used methods based on thymidine analogs, such as BrdU and EdU, in terms of accuracy and reliability.


INTRODUCTION
The neocortex constitutes the main part of the mammalian brain, and the location where the processing of all higher-order brain functions resides. Understanding its formation is one of the major interests in the field of developmental biology (Lodato and Arlotta, 2015). The neocortex develops from a stratified neuroepithelium, called the neural tube, into a complex structure of six horizontal layers of excitatory and inhibitory neurons (Matsuzaki and Shitamukai, 2015). Neurogenesis in the developing neocortex initiates when self-renewing neuroepithelial progenitors (NEPs) transform into apical and basal radial glia (RG) progenitor cells and start to produce neurons and intermediate neuronal precursors (Beattie and Hippenmeyer, 2017;Taverna et al., 2014). Since the discovery that RG constitute the progenitors of potentially all neurons in the vertebrate neocortex (Frederiksen and McKay, 1988;Hartfuss et al., 2001;Miyata et al., 2001;Noctor et al., 2004), a great effort has been focused in identifying their features and properties: how they coordinate in time and space to form the multiple layers of the neocortex; which signals control their fate; and how these signals orchestrate the correct balance between proliferation or differentiation during neurogenesis.
In principle, this balance can be robustly achieved via stochastic or deterministic cell decisions (Losick and Desplan, 2008). In brief, stochastic models assume a certain probability of differentiation that depends on the intracellular and extracellular signals that the cell is receiving. In this context, the fate at the single cell level is unpredictable and the balance between proliferation and differentiation is regulated at the level of the population (Teles et al., 2013). On the other hand, deterministic models of stem cell differentiation assume that the fate of the progeny is fixed and, therefore, the correct balance between the numbers of different types of neurons is achieved at the single cell level (Müller-Sieburg et al., 2002).
The dynamics of differentiation is often characterized based on the fate of the two daughter cells of a cell division relative to each other (Kosodo et al., 2004). This way, proliferating progenitors can perform pp ( progenitor-progenitor), pd ( progenitor-differentiated) and dd (differentiated-differentiated) divisions (Huttner and Kosodo, 2005). In this context, differentiation in the developing chick spinal cord (Míguez, 2015), in the zebrafish retina (He et al., 2012;Chen et al., 2012), epidermis (Clayton et al., 2007), airway epithelium (Teixeira et al., 2013), germline  and the intestine (Snippert et al., 2010) of mice follow a stochastic model. In these systems, progenitors can potentially perform each of the three types of division, and the corresponding rates are probabilistic and change overtime. On the other hand, the differentiation of RG in the mammalian brain has been shown to follow a deterministic asymmetric-only mode of division (Gao et al., 2014;Beattie and Hippenmeyer, 2017).
Several years ago, the group of Austin Smith showed that RG extracted from mouse developing neocortex can be successfully cultured in vitro (Conti et al., 2005). Driven by the multiple phenotypic similarities between neuronal precursors differentiated from embryonic stem cells in culture and RG, authors suggested that these neuronal precursors are the culture analogs to RG. In the same paper and driven by this observation, they also showed that in vitro cultures of RG could be established with fibroblast growth factor 2 (FGF2) and EGF as the key molecules that facilitate their expansion (Conti et al., 2005).
FGF2 is an extensively studied neurogenic factor for proliferation and differentiation of multipotent neural stem cells both during development and in the adult mouse brain (Kang and Hébert, 2015). FGF2 has been shown to be necessary for cell proliferation and neurogenesis in vivo, and to induce additional mitoses in progenitor cells in vitro (Raballo et al., 2000). In addition, stem cells from the adult mouse brain have been shown to proliferate and self-renew in vitro in the presence of FGF2 (Gritti et al., 1996). On the other hand, FGF2 stimulation have been shown to control the fate, migration and differentiation but not the proliferation of neuronal progenitors in vivo (Dono et al., 1998), whereas more recent studies do show an impact in promoting the cell cycle progression in cultures of rat glioblastoma cells (Baguma-Nibasheka et al., 2012).
From all these potential effects of FGF2, the specific features that facilitate the transition of RG from a non-expanding population in vivo that can perform only asymmetric pd divisions (and is, therefore, incompatible with progenitor cell expansion in number), to a self-renewing in vitro culture have not been quantitatively characterized in detail. In principle, this transition can be achieved by reducing the rate of neurogenesis, by promoting proliferative (at the expenses of asymmetric or symmetric differentiative) divisions, by increasing the proliferation rate (by shortening the cell cycle), by inducing cell cycle re-entry of quiescent progenitors (i.e. increasing the growth fraction), by reducing apoptosis (as a pro-survival signal), by inducing intermediate progenitors (that perform additional terminal divisions) or by shifting RG towards its less mature NEP phenotype (that perform pp divisions in vivo).
In this paper, we quantify the specific effects of FGF2 on key features of the proliferation and differentiation dynamics of RG that allow them to be cultured and expanded in vitro. To achieve this, we quantify values of cell numbers of RG and differentiated neurons extracted from mouse developing cerebral cortex and cultured in the presence of different FGF2 concentrations and at different time points. These values inform a theoretical framework based on a branching process formalism (Míguez, 2015) that provides average values of mode and rate of division of the RG population with temporal resolution. Our results show that FGF2 does not affect the rate of neurogenesis (i.e. the amount of differentiated neurons produced), it does not promote the NEP or intermediate progenitor phenotype and it does not affect significantly the apoptosis rate. On the other hand, FGF2 does promote symmetric pp divisions, it increases the growth fraction and shortens the average cell cycle length. These three key effects when combined, strongly facilitate the propagation and expansion of the culture.
In addition, discrepancies between predictions for the cell cycle length and growth fraction using several methods in our study pointed us to compare the accuracy of several common methodologies used to measure cell cycle features. To do that, we use a numerical model to show that methods based on cumulative thymidine analogs (such as EdU and BrdU) are not accurate in conditions of variable differentiation dynamics. On the other hand, the method based on branching process formalism performs better when mode and/or rate of division are changing, which is the case in our RG cultures and many in other in vivo developmental systems. In addition, the branching process method is superior due to its temporal resolution, robustness, minimal interference with cell homeostasis and simplicity of use.

FGF2 stimulation increases the growth rate of cultures of RG by shortening the length of the cell cycle
To initially test how the dynamics of growth and differentiation of RG in vitro is modulated by FGF2, cells derived from the developing neocortex of mouse embryos at E11-E11.5 were extracted, plated and cultured following standard protocols (Hilgenberg and Smith, 2007). Starting at 24 h post-plating (hpp), samples were fixed at three different time points and stained with Hoechst (Fig. 1A). Quantification of the number of cells in a field of view of fixed dimensions (0.6 mm×0.6 mm) using an automated segmentation tool developed in house (see Materials and Methods) is shown in Fig. 1B for two culture conditions, SC and SC+FGF, where the standard culture media is supplemented with an increased concentration of FGF2 ligand (see Materials and Methods). In both conditions, the number of cells increases, but the growth is only statistically significant (P<0.05) in SC+FGF conditions.
To study in detail how FGF2 affects the length of the cell cycle of the cycling progenitors, we performed 5-ethynyl-2′-deoxyuridine (EdU) cumulative labeling experiments to measure changes in the length of the average cell cycle. BrdU (Nowakowski et al., 1989), EdU (Salic and Mitchison, 2008;Buck et al., 2008) and other thymidine analogs constitute the most-used tools for estimating the cell cycle length of cells in many contexts (Alexiades and Cepko, 1996). The method is based on the replacement of endogenous thymidine during DNA synthesis with traceable compounds (Takahashi, 1966;Takahashi et al., 1996). The length of the average cell cycle is then inferred from the dynamics of the incorporation of these compounds into the DNA of cycling cells (Macdonald, 1970).
To estimate the average cell cycle length of the population, samples were cultured in the presence of EdU and then fixed at different time points (corresponding to different times of EdU incorporation). Combined nuclear Hoechst staining with EdU detection assay and immunostaining for Sox2 were used to identify all progenitors that had passed through S phase for each EdU incubation time.
The cell cycle length T and the growth fraction γ were calculated using the standard cumulative curve method based on linear regression (see Materials and Methods). Representative snapshots are shown in Fig. 1C-F. The resulting cumulative curves (Fig. 1E,F) reveal that γ remains at around 72% for both conditions tested, while T depends strongly on the culture conditions (T=35.2±3.2 h for SC, T=24.7±2.0 h for SC+FGF). In conclusion, our results show that FGF2 stimulation shortens the average cell cycle length in cultures of RG in vitro, while its effect in the growth fraction is not statistically significant. The previous section shows that FGF2 affects the rate of division.
To study the effect of FGF2 in the number of cells of each specific population of RG progenitors and differentiated neurons, we extracted the neocortex of mouse embryos at E11-E11.5 and plated cells at the same initial cell density in different wells. Next, cells were cultured under the two FGF2 conditions and samples were fixed every 2-4 h, starting at 24 h post-plating (hpp). Next, samples were stained using antibodies against Sox2 and Map2 to identify progenitors and differentiated cells, respectively. We then identified the fate of each cell based on the intensity of Sox2 and Map2 staining using our segmentation framework (see Materials and Methods).
Results are shown in Fig. 2A. Output provided by the segmentation script is plotted in Fig. 2B,C. Assuming the typical logistic growth model (Juarez et al., 2016) for proliferating cells in cultures, the corresponding sigmoidal curve fitting is also plotted (green, red and blue lines for RG, neurons and total cells, respectively). The data show that an initial regime of reduced change in cell numbers is followed by an increase in both cell types until the system reaches a regime where few new cells are being generated. Under both conditions, the amount of progenitors (green data points, green line) and differentiated cells (red data points, red line) increases with statistical significance (P<0.05) but the increase in progenitors is statistically more significant in conditions of SC+FGF (P=7.25E-09) that in SC conditions (P=7.60E-03).
In principle, this increase in the progenitor population could be explained by an increase in neuroepithelial progenitors (NEPs) (Beattie and Hippenmeyer, 2017;Taverna et al., 2014) or intermediate progenitors (Molyneaux et al., 2007), which emerge from asymmetric division of the RG and are able to perform a terminal dd division (Hutton and Pevny, 2011). Immunofluorescence against Tbr2, a marker for intermediate progenitors, shows no Tbr2-positive cells in the two culture conditions tested (Fig. S1C). This is in agreement with the effect of FGF2 in inhibiting the transition from RG to intermediate progenitor (Kang et al., 2009) (FGF2 is in the culture media in both experimental conditions: SC and SC+FGF). Quantification of immunofluorescence against Pax6, a well-characterized marker for RG (Suter et al., 2009) that is not present in NEPs (Elsen et al., 2018), shows that close to 100% of all Sox2-positive progenitors are also positive for Pax6 (Fig. S1A), suggesting that FGF2 stimulation does not promote the transition of RG back to neuroepithelial progenitors (Englund et al., 2005).
In conclusion, the increase in FGF2 concentration results in more RG and similar number of differentiated cells, showing that the population of cycling progenitors does not grow at the expense of the terminally differentiated cells.
Branching process formalism predicts variable mode of division that is affected by FGF2 stimulation The previous observation suggests that, apart from the changes in the cell cycle length, FGF2 may also be affecting the mode of division of the RG. It has been shown previously that the fate of differentiating RG can be modulated by FGF2, by changing the differentiation progeny of RG from neurons to glia (Qian et al., 1997). To quantify the effect of FGF2 in the mode of division, we took advantage of a branching process theoretical formalism developed by our lab (Míguez, 2015). In brief, the tool provides the average rates of each mode of division with temporal resolution simply based on numbers of progenitors and differentiated cells at different time points (see Materials and Methods).
Input data of the framework are the numbers of progenitors and differentiated cells, the rate of apoptosis and the growth fraction. To obtain the average rate of apoptosis, we performed immunostaining against anti-cleaved caspase 3 at three time points in the cultures under SC and SC+FGF conditions. Comparison between both conditions showed a very reduced rate of apoptosis that is not significantly affected by the addition of extra FGF2 (Fig. S1B).
Next, the apoptosis rate and the fitted values of cell numbers for progenitors and differentiated cells were used to estimate the average mode of division. This is provided by Eqn 1 in the Materials and Methods (Míguez, 2015) in terms of the difference between the rates of pp and dd divisions ( pp−dd). Results are shown in Fig. 3A. Interestingly, both conditions of SC (blue) and SC+FGF (red) show values of pp−dd≠0, which would correspond to the in vivo situation of asymmetric-only divisions pd=1 (as pp+pd+dd=1). In addition, the average rate of differentiation is not constant in time, with the maximum change in the differentiation dynamics occurring around 36-37 hpp, which is more prominent for SC conditions. Comparison between the two curves shows that the value of pp−dd predicted is higher when more FGF2 is present in the culture media, which corresponds with the higher increase in the number of progenitors observed in SC-FGF conditions (Fig. 2B,C).
To further validate the result that an increase in FGF2 increases the amount of pp divisions, we designed an experiment based on pulse and chase of EdU-labeled cells. To do this, we plated cells from mouse developing neocortex following the procedure explained in the Materials and Methods section. Next, cells were cultured under SC and SC+FGF conditions until 33 hpp. At this point, a 30 min pulse of EdU was applied to all samples. A number of samples were fixed at this time point (and labeled as 'Pulse' time point). The rest of the samples were washed with fresh culture media five times to remove the EdU (see Materials and Methods). These samples were cultured for another 15 h (corresponding to the predicted average T for SC+FGF conditions during this time, to ensure that labeled cells cannot cycle more than once in any of the culture conditions). Next, cells were fixed at this 'Chase' time point and stained using Hoechst, EdU and Sox2 immunostaining. Finally, the number of Sox2 + /EdU + cells at the time of the 'Pulse' (33 hpp) and 'Chase' (48 hpp) was quantified using our automated image analysis tool (see Materials and Methods). Results are shown in Fig. 3B,C. The number of progenitors labeled with EdU does not change significantly in SC conditions, consistent with a large proportion of asymmetric divisions (i.e. one EdU + RG produces two EdU + cells: one RG and one neuron; so the amount of EdU + RG remains constant) or a balanced ratio between pp and dd. On the other hand, in conditions of SC+FGF, we see a statistically significant (P<0.05) increase in the number of EdU + RG when comparing 'Pulse' and 'Chase' time points. This result shows that more RG originally labeled with the short EdU pulse divided and produced more RG when FGF2 was increased.
The length of the cell cycle is variable and shortens in response to FGF2 stimulation The branching process formalism also provides the average cell cycle length of the progenitors in the culture with temporal resolution (Eqn 2 in the Materials and Methods). This equation uses as additional input the value of the growth fraction γ, which can be indirectly obtained from the EdU experiments in Fig. 1E,F. A more direct estimation of the amount of quiescent progenitors, can be measured by immunofluorescence against KI67 at different time points (Fig. 4A) (Scholzen and Gerdes, 2000). The automated quantification of the number of Sox2 + cells that are also KI67 + in both culture conditions (Fig. 4B) shows statistically significant differences between SC and SC+FGF conditions, contrary to the results obtained with EdU cumulative curves in Fig. 1E,F. Under SC conditions, the growth fraction is around 55%, while the value in SC+FGF conditions is closer to 90%. This discrepancy between the EdU data ( Fig. 1E,F) and the KI67 immunofluorescence ( Fig. 4B) is discussed and studied in detail in the next section.
The value of the cell cycle length obtained as output of Eqn 2 is plotted in Fig. 4C, showing an average value of T that is not constant: a continuous decrease in cell cycle length is followed by an increase at later time points, and the minimum values for SC (around T=19 h) and SC+FGF (around T=10 h) conditions occur around 36-37 hpp. This reduction in the cell cycle length is in agreement with previous studies by the Dehay group showing that RG cells accelerate the cell cycle in response to FGF by shortening the length of the G1 phase (Lukaszewicz et al., 2002). In addition, these values are close to the values measured in vivo by Gao et al. (2014) and Beattie and Hippenmeyer (2017), who reported an average cell cycle length of 16-18 h in the temporal window corresponding to E11-E13.
In addition, we have observed that some cells in the culture organize in clusters, whereas others remain more isolated. This reorganization suggests that maybe there is heterogeneity in the differentiation dynamics (i.e. cells in clusters can retain apicobasal polarity and therefore behave differently from the rest of the cells in the culture). To test this hypothesis, we compared the change in EdU in the pulse-chase experiment (introduced in Fig. 3B) in cells in clusters versus non-cluster cells. Results are shown in Fig. S5, where the percentage of EdU + cells in both 'cluster' and 'non-cluster' locations, at both 'Pulse' and 'Chase' time points, and under both SC and SC+FGF conditions are plotted. Under SC conditions, the increase in EdU + cells in clusters and non-cluster cells is similar, suggesting a similar cell cycle length. Under SC+FGF conditions, the amount of EdU + cells is higher in the 'Chase' in cluster cells, but this difference may be due to an already higher number of EdU + cells at the 'Pulse' time point already in clusters. In fact, the difference between 'Pulse' and 'Chase' remains around 40%, suggesting that, also under SC+FGF conditions, the cell cycle is similar between clusters and non-cluster cells.
The branching process tool outperforms cumulative curve methods to monitor cell cycle dynamics Interestingly, and despite showing the same trend of shortening T with FGF2, the values of the cell cycle length predicted by the branching process formalism do not agree with the ones obtained by the EdU cumulative experiments in Fig. 1E,F. This discrepancy in cell cycle dynamics and in the growth fraction (Fig. 4B) led to us investigate the potential source of conflict between the cumulative method and the branching process tool by developing a numerical model of a generic differentiating stem cell population. A numerical analog of EdU is also simulated computationally, in such a way that cells in S phase are marked as labeled when EdU is present. The number of progenitors, differentiated and EdUpositive progenitors at each time point is used to calculate the average cell cycle length of the population using three widely used EdU-based methods: single cumulative curve (C1) (Nowakowski et al., 1989), dual cumulative (C2) (Shibui et al., 1989) and the pulse-chase (PC) method (Weber et al., 2014). The cell cycle is also calculated using the branching process (BP) method (Míguez, 2015) (Eqn 2 in the Materials and Methods). A detailed description of each method and how it is applied in this context is illustrated in Fig. S3 and explained in the Materials and Methods section. All predictions were then compared with the input value of T used for each simulation, to estimate the accuracy and reliability of each method.
The first scenario tested corresponds to homeostasis in the progenitor population ( pp−dd=0), constant value of T=20 h and no quiescent or apoptotic cells ðg ¼ 1; ; P ¼ 0Þ. These are the conditions defined by Nowakowski and co-workers when introducing originally the cumulative curve method (Nowakowski et al., 1989). Results of the analysis are plotted in Fig. 5A. Dots in Fig. 5B correspond to the prediction of the value of T for 10 independent simulations (crosses represent the average). We see that, for these particular settings, all four methods are able to predict the correct value of T (dashed line) within a 10% error margin, with both PC and BP performing slightly better than C1 and C2. Importantly, when comparing the individual values for the 10 simulations predicted by single and double cumulative curve methods (C1 and C2), there is a higher dispersion than in PC and BP methods. This means that a high number of repeats should be necessary to obtain an accurate value of T, and that the typical experimental design that involves only three independent repeats does not guarantee a correct estimation of the cell cycle. The same conclusions apply when considering growth of the population of progenitors, as in the case of RG reported here (pp−dd>0; Fig. S4A).
Variable cell cycle dynamics have been reported in many developmental systems (Míguez, 2015;Saade et al., 2013;Le Dréau et al., 2014;Takahashi et al., 1995;Calegari and Huttner, 2003;Calegari et al., 2005;Dehay and Kennedy, 2007;Mairet-Coello et al., 2012;Roccio et al., 2013;Arai et al., 2011;Iulianella et al., 2008;Locker et al., 2006). Fig. 5C shows the output of the numerical model when a variable value of T is used as input (with an average value T=20 h). Fig. 5D plots the quantification of the cell cycle in these conditions. In this situation, C1 predicts a much longer cell cycle that the average (49% error), while the C2 predicts a shorter cell cycle (24% error). Interestingly both PC and BP return a value much closer to the correct average, with less than 10% error. Again, the variability of the single cumulative C1 method (the one used in Fig. 1E,F and the most commonly used in the literature) is very high, making it unreliable when a small number of repeats are used (less than 10). Again, the same conclusions apply when considering conditions where the cell cycle changes while the population of progenitors is allowed to grow ( pp−dd>0; Fig. S4B).
The balance between differentiative and proliferative divisions has been shown to also change over time in many developmental systems (Saade et al., 2013;Míguez, 2015). We show here that, even in vitro, with cells growing in constant controlled conditions, the mode of divisions is highly non-constant (Fig. 3A). When we set a variable pp−dd in our simulations, we observe that, again, both single C1 and dual C2 cumulative methods fail and show high dispersion between independent samples (Fig. S4C). The same occurs when both mode and rate of division are allowed to change simultaneously (Fig. S4D). In these more realistic conditions closer to our experimental findings (variable mode and rate of division), the branching process equation predicts a value that is closer to the one used in the simulations, and the variability between samples is highly reduced.
In addition, our simulations show that the growth fraction measured with C1 method (the one used in Fig. 1) works well when the cell cycle and mode or division are constant, but it is not accurate when the parameters are variable. Therefore, the method based on Ki67 immunostaining (Fig. 4B) is more accurate under these conditions.
In conclusion, these results show that methods based on cumulative curve labeling are not suitable when proliferation and/ or differentiation rates are not constant. This, together with the reported effect of BrdU and analogs in lengthening the cell cycle (Levkoff et al., 2008), and the high dispersion when comparing sets of cells growing under the same exact conditions, could explain the discrepancy values of the cell cycle reported in Figs 1E,F and 4C. In addition, the error in the growth fraction measured in Fig. 1E,F versus 4A can be due to the same problems. Both pulse-chase (PC) and branching process (BP) perform well, while the branching formalism has the advantage of providing temporal resolution, as well as accurate values of the average mode of division during the experiment.
Values from the branching process analysis are able to reproduce the experimental data To test whether the values provided by the branching process formalism are correct, we took advantage of the same numerical model of the differentiating stem cell population introduced previously. Now, the model is informed with the values of initial number of cells as in the experiments (Fig. 2B,C), the values of T and pp−dd (Figs 3A and 4C), and the growth fraction γ and apoptosis measured in the previous sections ( Fig. 4B; Fig. S1B).
The model returns numbers of progenitors and differentiated cells that can be compared with the values of progenitors and differentiated cells measured experimentally. Accordingly, if simulations are able to reproduce the experimental observations, then the values of T and pp−dd obtained using the branching process formalism are correct. This comparison is shown in Fig. 6A,B, where the number of progenitors (thin blue lines) and differentiated cells (thin red lines) are plotted for 30 independent simulations using the same input values. Thick blue and red lines are the sigmoidal fittings of the experimental data. The good agreement between simulations and experiments under both conditions (SC and SC+FGF) suggests that the values of T and pp−dd predicted by the branching equations are correct.

DISCUSSION
A detailed analysis of the dynamics of vertebrate neurogenesis involves a careful characterization of the features that regulate the dynamics of proliferation and differentiation of RG during the generation of the mammalian cortex. One of its most striking features is the fact that RG are restricted to an asymmetric mode of division in vivo, as opposed to a more probabilistic scenario observed in other developmental systems (Saade et al., 2013;Fig. 5. The branching process tool outperforms cumulative curve methods. Cumulative curves and quantification for single cumulative (C1), dual cumulative (C2), pulse-chase (PC) and branching process (BP) methods for 10 independent runs of the numerical model for conditions of (A,B) constant and (C,D) variable cell cycle length. Each color corresponds to the same simulation analyzed using each framework (see text). The cell cycle is also calculated using the branching process (BP). Dots correspond to single runs of the model, crosses show the average value for the 10 independent simulations.
Míguez, 2015; He et al., 2012;Chen et al., 2012;Clayton et al., 2007;Teixeira et al., 2013;Klein et al., 2010;Snippert et al., 2010). FGF2 has been shown to facilitate the expansion of RG in vitro cultures, but the details of this process have not been studied. Our quantitative characterization of the role of FGF2 shows simultaneous effects in the growth fraction (Fig. 4B), in the mode of division (Fig. 3A) and in the length of the cell cycle (Fig. 4C).
The overall influence of each of these effects in the expansion potential of the RG culture can be assessed using our numerical model. To do this, we informed the simulations with the experimental values for SC, and quantified the increase in the number of cycling progenitors (as a measure of the potential of the culture to expand in size) after 22 h. Next, we substituted each of the predictions for cell cycle length, growth fraction and differentiation rate predicted for the SC+FGF2 conditions, individually or in combination. The increase in cycling progenitors for 30 independent numerical simulations for each condition is shown in Fig. 6E. Surprisingly, the analysis suggests that the most influential feature is not the differentiation rate or the growth fraction, but the change in cell cycle length. The changes in growth fraction or the mode of division do not significantly impact the culture in terms of cycling progenitors (by 1% and 9%, respectively). On the other hand, the effect of mode of division combined with the cell cycle increases the expansion rate (51%), when compared with the effect of the cell cycle alone (34%).
Several authors propose that the mode of division depends on the distribution of cell fate determinants during mitosis, the orientation of the spindle or the inheritance of the primary cilium or the different centrosomes (Taverna et al., 2014). It is possible that the apical-basal polarized structure of the RG, or their organization and orientation of the radial processes along the stratified neuroepithelium results in asymmetric inheritance of these cell fate regulators (Taverna et al., 2014). The loss of these polarizing features provided by the niche when cells are cultured in vitro may result in a probabilistic scenario where the fate of the two daughter cells is independent of each other and all three modes of divisions are possible, similar to neuronal progenitor cells and other developmental systems (Saade et al., 2013;Míguez, 2015;He et al., 2012;Chen et al., 2012;Clayton et al., 2007;Teixeira et al., 2013;Klein et al., 2010;Snippert et al., 2010). In fact, early studies in the mouse neocortex suggest that the model that best fits the clone distribution assumes that the fates of the daughter cells are independent of each other (Cai et al., 2002). In this situation, the branching process framework is able to estimate the rates of each of the three modes of division (Míguez, 2015). This prediction for the case of RG in culture is shown in Fig. 6C,D, where we can see that the predominant mode of division is pp (green). This symmetric mode of division is even more probable under conditions of SC+FGF, at the expense of reductions in pd and dd.
A detailed analysis of the dynamics of vertebrate neurogenesis involves a careful characterization of the rate of division. The mostdirect method to measure the cell cycle length requires the monitoring of the time between consecutive mitotic events at single cell resolution (Sigal et al., 2006;Wilcock et al., 2007). Unfortunately, owing to the high degree of variability, many cells in a population need to be sampled, segmented and tracked simultaneously to obtain an accurate value, even when dealing with clonal samples (Sandler et al., 2015). Therefore, the most-used approach is that of thymidine analogs, but this has several drawbacks: it can be toxic and mutagenic (Duque and Gorfinkiel, 2016), and can affect the normal dynamics of cell proliferation (Levkoff et al., 2008) by lengthening the cell cycle. In addition, choosing the correct mathematical analysis and interpretation of the experimental data is not straightforward (Johansson et al., 1999). Authors have proposed several approaches, such as linear fitting (Begg et al., 1985;Høyer et al., 1994), nonlinear fitting (Johansson et al., 1994;Weber et al., 2014), or the use of deterministic (Lee and Perelson, 2008) or stochastic models (Zilman et al., 2010). Depending on the method used, the same input data results in very different predictions for the average duration of the cell cycle (Ritter et al., 1992). Owing to these limitations, BrdU and analogs have been referred to as 'one of the most misused techniques in neuroscience' (Taupin, 2007).
Our results show that methods based on cumulative incorporation of thymidine analogs perform well in conditions of constant proliferation and differentiation, but they are not designed to study systems where the cell cycle changes over time, which is potentially the case in many developmental systems. Under these conditions, the branching process formalism and the pulse-chase outperform cumulative curve methods. On the other hand, the pulse-chase method requires experiments that are longer than the cell cycle length, so an estimation of the value of the cell cycle has to be known beforehand. In addition, the well-known toxic effect of the labeling agent for such long periods of time may affect strongly the normal cell cycle progression, probably by enlarging its real value (Levkoff et al., 2008;Duque and Gorfinkiel, 2016). A clear advantage of the branching process is that it does not involve manipulation of the samples before fixation, so there is no interference with the normal progression of the cell cycle. In addition, the branching process formalism also provides the correct value of T with temporal resolution and the measurement of the average differentiation rate (also with temporal resolution). Some of the more important limitations of the branching process tool are that it requires absolute values of numbers of cells and not relative amounts (thus, data from flow cytometry is not useful). In addition, it can only be used in situations where differentiated cells are terminally differentiated.
Several studies have shown that the length of G1 phase increases progressively when neurogenesis starts, resulting in a overall increase in cell cycle length (Takahashi et al., 1995;Calegari and Huttner, 2003;Calegari et al., 2005;Dehay and Kennedy, 2007;Mairet-Coello et al., 2012;Roccio et al., 2013). Alternatively, other studies show that the cell cycle length is shorter in neurogenic divisions, compared with proliferative divisions (Arai et al., 2011;Saade et al., 2013;Le Dréau et al., 2014;Iulianella et al., 2008;Locker et al., 2006), owing mainly to a shorter S phase. The potential effects of FGF2 in each phase of the cell cycle has been studied previously (Hitomi and Stacey, 1999). Our results show that FGF promotes pp divisions and shortens the cell cycle, consistent with the findings of Lukaszewicz et al. (2002), and with the hypothesis that proliferative divisions have a shorter cell cycle, maybe due to a shorter G1 phase, similar to the effect of insulin-like growth factor (Mairet-Coello et al., 2009;Hodge et al., 2004).
Based on this, two of the consequences of FGF stimulation (shortening T and inducing pp divisions) may correspond to two aspects of the same role of FGF. For example, FGF induces pp divisions that have been shown to be faster than neurogenic divisions (Calegari et al., 2005). On the other hand, the lengthening of the neuroepithelial cell cycle has been shown to be sufficient to trigger differentiation. Based on this 'cell cycle length hypothesis' (Calegari and Huttner, 2003), the shortening of the cell cycle mediated by FGF may result in more pp divisions. The link between cell cycle length and mode of division is a very interesting topic, and our branching process tool is a powerful tool with which to study this interplay (Míguez, 2015).
Another potential explanation of our observations takes into account the mutual inhibitor link between FGF and retinoid signaling (del Corral and Storey, 2004). As retinoid signaling promotes the expression of inhibitors of CDKs, as well as neurogenic genes, its repression by increased levels of FGF can explain the triple effect of FGF in promoting proliferation, accelerating the cell cycle and even increasing the growth fraction.
An additional physiological interpretation of our observations involves the reported ventralizing potential of FGF (Pilz et al., 2013). It has been shown that FGF drives expression of proteins and transcription factors typical of a ventral radial glial identity, such as progenitor cells in the developing ganglionic eminence.
Here, a large number of sub-apical progenitors show selfrenewing capabilities, and their cell cycle shortens with each generation. The fact that FGF results in cell cycle acceleration and more self-renewing divisions is consistent with the ventralizing potential of FGF.
In addition, previous studies in chick spinal cord have shown that FGF has an heterogeneous effect, inducing incomplete interkinetic nuclear migration and fast non-apical divisions in a subset of progenitors. Unfortunately, our approach based on fixed samples and immunofluorescence does not allow us to monitor individual cellular responses to investigate potential heterogeneous responses in the pool of RG.
Another interesting approach is to study the interplay between the two key signals (FGF and EGF) that have been shown to be key to establish the self-renewing expanding RG culture (Conti et al., 2005). A currently ongoing project of our lab is to take advantage of the framework used here (in vitro culture, automated image analysis, modeling and branching processes) to establish how FGF and EGF interplay to fully inhibit differentiation and promote self-renewing divisions in tissue cultures of RG (Conti et al., 2005).

Conclusions
The culture and differentiation of RG cells in vitro provides a very good framework to study basic features that orchestrate the formation of the mammalian neocortex. In brief, the system provides a well-controlled environment where the effect of signaling molecules and other conditions can be tested reliably, while providing easier manipulation and imaging compared with studies performed in vivo. We use this framework to study the features that promote the expansion of RG in culture driven by FGF2. Our combined experimental/computational/theoretical approach can be also used to test the effect of other signaling networks by quantifying the cell cycle and mode of division after ligand stimulation or small molecule inhibition, after a comparison with a control culture.

MATERIALS AND METHODS
Preparation and culture of dissociated mouse cortical RG Cells were obtained from mouse embryos of the C57 BL/6JRCC line at E11/ E11.5, following standard methods described previously (see Hilgenberg and Smith, 2007). The initial time point is labeled as 0 h post-plating (hpp) and it is used as the reference point for our experiments. Briefly, after careful removing of the meninges, the cortex was isolated and placed in Hank's Buffered Salt Solution free of Ca 2+ and Mg 2+ (HBSS, ThermoFisher 14185). Next, samples were mechanically disaggregated using Pasteur pipettes and plated in coverslips treated with nitric acid and fibronectin at 10 µg/ml (Fisher Scientific; 15602707) to facilitate cell adhesion. Cells were plated at constant density (250,000 cells in each M24 well) for all experiments in Neurobasal medium without L-glutamine (ThermoFisher 21103-049), Glutamax (ThermoFisher 35050-038), B-27 (ThermoFisher 17504-044), penicillin, streptomycin and anti-mycotic (concentrations standard for cell culture). Media were complemented with 0.02 ng/µl of recombinant murine EGF (PeproTech 315-09, lot number 0517179-1) and 0.02 ng/µl of human basic FGF (PeproTech 100-18B, lot number 0311706-1). This culture medium is referred to as the standard culture (SC) condition in our study. Cells are allowed to rest for 1 day in the incubator to recover from the dissection process. At 24 hpp, the culture medium was changed for fresh SC media or for SC media complemented with additional human basic FGF at a final concentration of 0.06 ng/µl. This culture condition is labeled as SC+FGF in this study. All experimental protocols were carried out in accordance with the guidelines of the European Communities Directive (2012/63UE) and Spanish legislation (RD 53/2013).

Statistical and data analysis
One-way ANOVA test was used to measure statistical significance between different time points. Cell cycle values in Fig. 1E,F were obtained after linear regression of the four first data points. Rates of quiescence in Fig. 4B were obtained from the mean value of the four last points. Slope error was calculated using a linear fitting with values of the average plus standard error and another linear fitting with values of the average minus standard error to obtain the difference in the slope between these two values. Quiescence error was the standard error of the four last points, and the T error was derived from the error propagation of the previous values. Three-parameter sigmoidal fitting was used to fit data from Fig. 2B,C. Sample size for all experiments was at least four. Unless specified, error bars represent the standard error of the mean, calculated using the error propagation. Confidence intervals for the prediction of T and pp−dd in Figs 3A and 4C were calculated by taking half of the value of the maximum error for P and D in Fig. 2B,C. Next, we fitted these data points to three-parameter sigmoidal functions. Finally, these curves were used as inputs to the branching process equations to obtain the 50% confidence intervals plotted. All curve fitting and statistical analysis were performed using Matlab (The Mathworks) and Julia programming language (Statistics package).

EdU cumulative curve
A cumulative curve of the thymidine analog 5-ethynyl-2′-deoxyuridine (EdU) incorporation was created using a Click-iT Plus EdU Alexa Fluor 647 Imaging Kit (ThermoFisher, C10640). Briefly, EdU was added around 24 hpp at 2 μM. Cells were then fixed at increasing times of EdU exposition. Staining of EdU-positive cells was performed based on previously published protocols (Harrison et al., 2018). Immunostaining for Sox2 was used as a standard marker for RG progenitors (Beattie and Hippenmeyer, 2017). Later, the number of cells positive for both Sox2 and EdU was quantified using automated image processing. To calculate the cell cycle length, the percentage of progenitor cells that have incorporated EdU was plotted against the hours of EdU incorporation. The saturation value at long incubation times was used to calculate the growth fraction γ. This value was then used to calculate the average cell cycle using linear regression at short EdU accumulation times (see Fig. 1).
EdU pulse-and-chase experiments Cells were exposed to a short pulse of 30 min of EdU at 36 hpp. 'Pulse' points were fixed at this time point. 'Chase' points were washed three times with fresh medium and were fixed 15 h after the 'Pulse' time point. The number of EdU-positive/Sox2-positive cells is quantified for both 'Pulse' and 'Chase' time points for both conditions using automated image processing.

Branching process formalism
Our lab has developed a method to measure the dynamics of proliferation and differentiation that does not depend on cumulative labeling with thymidine. Instead, it uses a branching process formalism, where cell cycle, growth fraction and mode of division are defined as independent variables, to derive analytical equations that provide the average values of proliferation and differentiation of the population based only on the numbers of proliferative, differentiated, quiescent and apoptotic cells at different time points. A scheme of the method is shown in Fig. S3D, and an example of its experimental implementation can be found in Míguez (2015). In brief, the branching process tool works with a minimum of two time points, and it does not require data points to be equispaced in time. In terms of repeats, depending on the interests of the works, it would be more useful to include more repeats to obtain more-reliable values, or more data points to obtain a more-detailed time evolution of the values. The number of data points needed to obtain a reliable temporal evolution depends on how fast cell numbers change (as a rule of thumb, we suggest a minimum of five time points each time the number of total cells doubles). In brief, if the user has a series of measurements of P and D (n experiments) for different time points (t 0 , t 1 , …, t m ), they can simply use directly the mean of these data points as input of the equations to obtain discrete values of proliferation and differentiation rates between the time points. Another option is to fit the experimental values (using the mean and the standard deviation) to a representative polynomial curve (as we do here) and obtain a more smooth time evolution of the proliferation and differentiation rates.
To obtain these values, samples are allowed to develop without interfering with the normal dynamics of the cells, and then are fixed at different developmental times. After fixation, the amount of cells in each state is quantified by antibody staining to distinguish progenitors (P) (Graham et al., 2003), differentiated (D) (Míguez, 2013) and the number of progenitors undergoing apoptosis (Ø P ) (Blanchard et al., 2010). The growth fraction γ is obtained using double immuno-labeling against Sox2 and Ki67 (Scholzen and Gerdes, 2000).
These values (quantified using the automated quantification described) are then fitted to sigmoidal functions (green and red curves in Fig. 2B,C) that are used as inputs of the following equations for the mode and rate of division, which correspond to a generalization of equations presented by Míguez (2015) updated to account for a potential reduction of the progenitor pool: pp À dd ¼ 1 þ P ððDD=DPÞ À 1Þ ðDD=DPÞ þ 1 (Eqn 1) where pp and dd correspond to the rate of symmetric proliferative and differentiative divisions, respectively. ΔP=P t −P 0 and ΔD=D t −D 0 correspond to the number of progenitors and differentiated cells generated in a given window of time Δt=t−t 0 . The value pp−dd goes from 1 (all divisions are symmetric proliferative) to −1 (all divisions are symmetric differentiative). The value, pp−dd=0 corresponds to maintenance of the progenitor pool, either via asymmetric pd divisions or via a balance between symmetric proliferative and differentiative divisions (the model cannot distinguish between these two scenarios, as they are mathematically equivalent). ∅ P is the rate of cell death of the progenitors pool, obtained using double immunolabeling against Sox2 and cleaved caspase 3 (see Fig. S1B). This reduced value of apoptosis rate (assuming that most cell death occurs via apoptosis) is consistent with estimations from in vivo experiments (Cai et al., 2002).

Image acquisition and analysis
Samples were imaged using a confocal microscope (AR1+) that has a high speed of acquisition and sensibility coupled to an inverted microscope (Eclipse Ti-E, Nikon) with a 20× objective and a resolution of 1024×1024 pixels. The field of view was set to 0.6 mm×0.6 mm. In brief, image processing and analysis ( performed in Fiji; Schindelin et al., 2012) was based on the segmentation of nuclei and the classification of each cell as progenitor, differentiated neuron, quiescent or apoptotic based on the intensity of the fluorescence staining of each marker. A large number of cells (around 10 5 cells) was processed for each data point to minimize the effect of variability and heterogeneity of the samples. The sequence of processing algorithms and filters is as follows. 1. Definition of the Kernel Radius (KR) that sets the size of the region used for calculations and filter processing. Several KR sizes were tested (values from 1 to 5 pixels). The final KR was fixed as 2.5. 2. A local thresholding is applied to remove background based on the median intensity as cutoff value (radius=8×KR). 3. To remove breaks and holes inside the objects generated by the previous filter, the following sequence of filters is applied to enhance the definition of the boundaries of each object: Gaussian Blur filter, Maximum Filter, Median filter and Unsharp Mask filter (radius=KR). 4. The resulting image is binarized using the median value as threshold. 5. Euclidean distance mapping (EDT) is performed in the binary image to generate seeds that are used by a flood-fill algorithm to define the boundaries of each object (Kang et al., 2010). 6. Finally, all objects are fitted to ellipses for posterior analysis. Ellipses smaller than 4×π×KR 2 are discarded from the analysis. The specific features of each staining requires a different set of processing filters to enhance signal for each channel.
1. Map2: double sequential thresholding to extract foreground information (cutoff 1=mean, cutoff 2=median); morphological opening to remove neurons fibers (structuring element: lines at different angle with a length of 2×KR); Gaussian filter to remove noise (radius=KR). 2. Sox2: double sequential thresholding to extract foreground information (cutoff 1=mean, cutoff 2=median); morphological opening to select only nuclei with minimal size (structuring element: circumference of radius equal to 2×KR); Gaussian filter to remove noise (radius=KR). 3. EdU: single thresholding to extract foreground information (cutoff=median); morphological opening (structuring element: circumference of radius equal to 2×KR); Gaussian filter to remove noise (radius=KR). 4. Pax6: single thresholding to extract foreground information (cutoff=mean); morphological opening (structuring element: circumference of radius equal to 2×KR); Gaussian filter to remove noise (radius=KR). 5. Cleaved caspase 3: double sequential thresholding to extract foreground information (cutoff 1=mean, cutoff 2=mean+plus s.d.); morphological opening to select only nuclei with minimal size (structuring element: circumference of radius equal to 2×KR); Gaussian filter to remove noise (radius=KR). 6. Ki67: single thresholding to extract foreground information (cutoff=mean); morphological opening (structuring element: circumference of radius equal to 2×KR); Gaussian filter to remove noise (radius=KR). Finally, the identity of each ellipse was established based on the number of pixels above threshold in each channel. For the Map2, this area was set to at least 15%; for the rest it was set to 1%. A subset of cells was both Sox2− and Map2−, and had a nucleus that was much larger that Sox2+ or Map2+. As these were not RG or differentiated neurons, they were not taken into account in the study.

Numerical simulations of cell populations
We developed an in silico phenomenological numerical model of a generic differentiating stem cell population that simulates cycling progenitors that can either proliferate, differentiate, enter quiescence or enter apoptosis based on rates and probabilities provided by the user. Each cell has the following features: length of its cell cycle (T), current phase of cell cycle, time since birth (age) and fate ( progenitor, quiescent, differentiated or apoptotic). Values of cell cycle length, mode of division, quiescence and death rate can be kept constant throughout the simulation or can be set to change each timestep. For each cell in the simulation, the cell cycle length, the probability of differentiation at the end of the cell cycle, the probability of entering quiescence and the probability of entering apoptosis are obtained from gamma distributed values with a mean value set by the input parameters (T, pp−dd, γ, and ∅ P ) and standard deviation of 30% of the mean to mimic intrinsic cell-to-cell variability and intrinsic noise in a clonal population (León et al., 2004) (other values of the standard deviation from 10% to 50% provide similar results).
A scheme of how the population is defined and develops over time is shown in Fig. S2. Parameters of the simulation are: the number of initial cells m, the average cell cycle T at each time point (defined as T ¼ P T i=n, n being the number of cells at time t), the fraction of cycling progenitors (or growth fraction) γ, the rate of apoptosis of progenitors ∅ P and the length of the experiment t end . The age of each cell is defined as the time since its birth, and the type corresponds to its characteristic as progenitors (P, cycling cells), differentiated (D, non cycling cells), quiescent (Q, non cycling progenitors) and apoptotic (dying cells). The simulation takes palace as follows: an initial set of unsynchronized progenitor cells are allowed to cycle following the different phases of the cell cycle: from G 1 to S G 2 finally M phase. Upon division, the two resulting daughter cells either remain as progenitors ( pp division), become terminally differentiated cells and stop cycling (dd division), or one remains as a progenitor while the other differentiates ( pd division). For simplicity, the cell cycle is divided into just three main steps of equal length: G 1 , followed by S and finally followed by G 2 +M (T=T G1 +T S +T G2M ). Changes in the cell cycle length affect all phases of the cell cycle identically (simulations where the phases are of different length and where changes affecting different phases of the cell cycle in different ways show equivalent results).

Simulations of cell cycle determination methods
The previous model is then adapted to perform a computational analog of one or two thymidine compounds. At any time in the simulation, EdU can be added to the cells, so cells undergoing S phase will be labeled as 'positive' and will remain as positive throughout the rest of the simulation. The input parameters of the model are varied to simulate different dynamics of a population of cells in different conditions, in terms of quiescence, apoptosis, cell cycle length and differentiation rate. For each condition tested, we perform four measurements of the cell cycle based on the following methodologies:

Cumulative curve method
This technique has been extensively used both in in vitro and in vivo situations to quantify the rate of cells in the population entering S phase (Martinez-Morales et al., 2010;Le Dréau et al., 2014). A scheme of the method is shown in Fig. S3A. In brief, a nucleoside analog is added to several identical samples that are fixed and stained at different times. Labeled cells in all samples are quantified using microscopy or flow cytometry. The ratio of progenitor cells that are labeled for each sample is plotted, and the values corresponding to the cell cycle length T are obtained from the slope of a linear regression fitting of the data at short exposure times. In addition, the fraction of cycling progenitor cells γ, or growth fraction, can be estimated from the rate of labeled cells after long exposure times. This method, when combined with dyes to measure DNA content, can be used to determine the length of the different phases of the cell cycle (Dolbeare and Selden, 1994).

Dual cumulative curve method
This method combines dual staining with thymidine analogs (Salic and Mitchison, 2008). It also provides the possibility of fixing all samples simultaneously to ensure that quantification is always performed at the same developmental time. In addition, it can also provide some positional information of regions in a given tissue where cells cycle at different rates (Shibui et al., 1989;Bradford and Clarke, 2011). However, it requires a more-complex experimental design, and it may also result in increased toxicity. In addition, it does not provide information about the growth fraction. The method (Fig. S3B) involves an initial labeling agent administered to all samples simultaneously, and a second agent administered at different time points. All samples are collected at the same time, and they are stained for both labeling agents. The amount of cells that are double positive over time for the two different thymine analogs is plotted, and the average length of T and T S can be obtained using linear or nonlinear regression (some corrections regarding the potential differential incorporation of both agents are required).

Pulse-chase method
Both previous methods rely on long-term exposure of the samples to nucleoside analogs, which can result in toxicity effects. Alternative, a short pulse can also be applied (Weber et al., 2014) to label only cells that were in S phase at a given time. The population of positive cells is then 'chased' in the different samples by fixing and staining at different times. Several variations of this method have been developed. A commonly used technique is to stain cells in mitosis (using immunofluorescence against phospho-histone 3), or using a second thymine analog in S phase to chase cells that have re-entered in a new S phase. A scheme of the method is shown in Fig. S3C.
The ratio of double-positive cells in the different samples is plotted over time, and the average value of T corresponds to the time between the pulse and the maximum number of double-positive cells in the population. The slope of the curve at shorter time scales can be used to calculate the length of S phase. Measurements of the cell cycle using this method require significantly longer experiments than the two previous methods.

Branching process method
The number of cells and their fate as progenitor, differentiated, quiescent or apoptotic cells is recorded at each time point during the simulation. These values are then used as input of the branching process (Eqn 2). The average value is then plotted for each condition tested.