Negative Autoregulation by FAS Mediates Robust Fetal Erythropoiesis

Tissue development is regulated by signaling networks that control developmental rate and determine ultimate tissue mass. Here we present a novel computational algorithm used to identify regulatory feedback and feedforward interactions between progenitors in developing erythroid tissue. The algorithm makes use of dynamic measurements of red cell progenitors between embryonic days 12 and 15 in the mouse. It selects for intercellular interactions that reproduce the erythroid developmental process and endow it with robustness to external perturbations. This analysis predicts that negative autoregulatory interactions arise between early erythroblasts of similar maturation stage. By studying embryos mutant for the death receptor FAS, or for its ligand, FASL, and by measuring the rate of FAS-mediated apoptosis in vivo, we show that FAS and FASL are pivotal negative regulators of fetal erythropoiesis, in the manner predicted by the computational model. We suggest that apoptosis in erythroid development mediates robust homeostasis regulating the number of red blood cells reaching maturity.

initially gating the cells at low Ter119 levels and then finding minima separating two peaks on the histograms of CD71 expression levels (these histograms along with the CD71 expression level histograms for all expression level of Ter119 are shown for all time points in Supplementary Figs. 1-14). By definition adopted here, this line (horizontal) separated states 1 and 2. The second line (vertical) was drawn similarly, by assaying Ter119 expression at high CD71 levels (not shown). This line separated states 2 and 3. The final line (horizontal) designed to separate states 3 and 4, was harder to identify, since no obvious separation between low and high levels of CD71 expression was seen for high Ter119 expression levels. However, two self-consistent observations indicated that there is indeed a difference between (high Ter119, high CD71) and (high Ter119, low CD71) stages. First, by looking at the data for embryonic days 13 through 14.5, one can see that there is a relative slowing or 'lingering' of the erythropoietic wave front at high CD71 and high Ter119 levels. Additionally, and more significantly, there is an identifiable separation on the erythropoietic wave, after which cells cease to decrease their size and change other characteristics ( see Sup. Fig. 15), indicating mature state representing state 4 in our analysis. One can draw a line outlining the position of the front during the slow-down of the wave or another line indicating where cell size change becomes undetectable. Interestingly, these two lines are localized to almost the same region, leading to a more defined separation between (high Ter119, high CD71) and (high Ter119, low CD71) states, and thus allowing us to introduce separation between states 3 and 4 ( Supplementary Fig. 15). We have thus drawn a second horizontal line separating these states. All three separation lines described above are shown in each of the panels in Supplementary Figs. 1-14. Some changes in absolute positions of the lines should be attributed to somewhat different flow cytometry settings adopted for different experiments corresponding to different time points. The next step was to weave together a variety of data points which were taken at various stages of the erythropoietic differentiation. One should remember that the erythropoieitc wave is a result of amalgamation of data from different embryos sacrificed at different estimated periods of their developmental time. We assigned approximate times to each of these data sets, placing them in evenly spaced intervals (0.5 embryonic days), and then later tested the sensitivity of the results to this assumption (Section 4). The time course thus represented an overall dynamic profile of erythroid cell differentiation, as defined by the acquisition and loss of cell surface markers CD71 and Ter119, with correlation to cell size.   Figure 15. Determination of the boundary between States 3 and 4 E14.5 fetal liver cells were labeled with antibodies against CD71 and Ter119. Viable cells (negative for the DNA dye 7-AAD) were selected (A). Ter119-high cells were gated and divided into regions of decreasing CD71 expression labeled 'a' to 'f' and color-coded as indicated (B). Histograms of Ter119 expression (C), side-scatter (SSC, D), and forward scatter (FSC, E) were constructed for cells in each region. These indicate a gradual increase in the mean Ter119 expression, reaching a maximum in regions 'e' (light blue) and 'f' (violet). Regions 'e' and 'f' also have similar histogram characteristics with respect to SSC, a function of cellular complexity, where they both reach a minimum. The FSC histograms demonstrate two populations, related to the cell-cycle status of the cells. Cells in S-phase and G2/M are likely to make up the higher FSC population. The position of the peak in the higher FSC population is related to cell size, and is the lowest for cells in regions 'e' and 'f'. Cells in regions 'e' and 'f' also have the highest proportion of cells in the lower FSC population. Taken together, cells in regions 'e' and 'f' represent the most mature erythroid precursor cells with relatively uniform characteristics. We therefore chose cells in regions 'e' and 'f' to constitue state 4, and the boundary between regions 'e' and d' in this figure to represent the boundary between states 3 and 4.
Each state now represents a proportion of the overall cell population, determined by the gating coordinates of Ter119 and CD71 expression (Table 1). Due to rounding inherent to the analysis in the WinMDI software, the data was re-normalized, making a small adjustment on the magnitude of each of the states. The results are shown together in Table 1, which was used for the subsequent analysis. We have also altered these data somewhat to test the sensitivity of results to various alterations of this dataset, as explained below. Based on these experimental data, we wanted to determine the most likely underlying mechanisms of interaction between cells in different states. We approached this problem by testing in a comprehensive fashion different hypotheses, each assuming different number and types of interactions between cells in states 1 through 4. These hypotheses were expressed as mathematical models in the form of ordinary differential equations (ODEs). We have thus explored a wide space of potential cell-cell interaction 'topologies', limited only by the combinatorial possibilities. As explained more in detail below, these models were tested for their ability to both represent the experimental data accurately, as well as remain resistant to uncertainty in model parameter values. The main assumption of our analysis was that interactions between cells in different states frequently present in 'good' models ('goodness' of a model will be further defined below).are also present with high probability in the real system. Furthermore, to deal with uncertainty in the experimental measurements, we identified the cell-cell interactions that are robust to variation of times assigned to specific measurements, and to other alterations to the basic experimental data set. The overall analysis was manifested as a large scale computational process that iteratively constructed the hypothetical models, and tested them for suitability for representing the true control architecture of red blood cell differentiation. The vector field for every system of equations is thus a balanced input-output relationship governed by the nonlinear rates of transition, hi,i+1. The system is mass-balanced, i.e. we assume that the total number cells at each time point remains unity. We assume that the only driving force is the initial conditions, as there are no exogenous inputs nor is there exit from the system, i.e., the cells in state 4 are assumed to be in the final differentiation state on the time scale considered.. Further constraints were added in subsequent analysis to retain non-negativity of both states and transitions, as described below.  always non-negative. The transition function hi,i+1 is also influenced by a second type of interaction, due to feedback or feedforward regulation of cells in the source state, xi , by cells in any other state, xj . The effect of such an interaction on the transition function may be either positive or negative, and is proportional to the product of the cells in the two interacting states, xixj. In the special case where j=i, the interaction is autoregulatory and is proportional to xi 2 . Twelve parameters, pk (k ∈ 4-15), represent the rate constants for each of these potential interactions, which are classified as feedforward [ff(j, Supplementary Table 3. The transition function hi,i+1 is therefore the sum of all these potential interactions (Eq. 2 below). Each of the feedforward or feedback inter-state interaction parameters (p4-p15), can be positive or negative, to be determined by regression for each model topology, and thus are of sign that can vary with the model topology. Together, these parameters contribute to each overall state transition, and define the "connectivity" or "structure" of each hypothetical model.

Equation 2. Generalized Transition Function
, 1 ( , ) The level of abstraction for the generalized transition function is partially based on the fact that the experimental dynamics are simple, and should not require long series and sums in order to reproduce them computationally.

Generation of different 'model topologies'
In all, in our analysis, 298 different possible differentiation network models and the corresponding vector fields were constructed by 'turning on/off' each of 12 different parameters regulating up to three nonlinear state dependencies (representing cell-cell interactions between cells in states 1-4) that we assumed could exist simultaneously in one model.
The assumption that the number of regulatory interactions is limited to three is based both on the desire to limit the combinatorial explosion in the number of possible state dependencies, and on the more general argument that multiple and dense interactions between elements of any network can lead to instabilities (2). In a similar fashion, a recent network reconstruction analysis considerably limited the number of possible links in potential networks, obtaining nevertheless excellent results (3). More precisely, the model networks were constructed and investigated as follows. As there are 4 distinct erythropoieic states (as defined above), there can exist 12 different potential interactions between these states (each state can interact with 3 others). These interactions can be both positive and negative, further diversifying possible network topologies. We opted to allow the sign of each parameter to be determined by data regression, and allowed only one, two, or three total connections to exist simultaneously in any individual model. This reduces our total network model number (without explicit specification of the sign of each interaction) to:

Model fitting to the experimental data
Each vector field varying in the structure of its transition functions was generally described by templates of 1's and 0's which signal to the algorithm to either vary the parameters to minimize an objective function (Equation 4) consisting of similarity and penalty functions, or remain zero, thereby defining a characteristic regulatory scheme. The similarity function ( , ) is the squared difference between the experimental data and each fitted trajectory. It represents the sum of the squared difference between the experimental data and the fitted curves for all 4 states ( f is a single value -the sum of the squared residuals for all four curves at once, for a total of 28 points). This fitting was performed for each separate model. The fitting was performed using different methods, including Gradient Descent and Random Parameter Search with similar results, but the more efficient algorithm that was also more effective in terms of avoiding local minima in optimization was Direct Version contains the optimum, and then divides this region in the next iteration. This process continues until the parameter space is optimally defined, i.e. the function optimized approaches a constant minimum in a small region of parameter space. This algorithm was used as written, with only minor modifications as to how parameters are sent to the module. In addition to minimizing the difference between the experimental data and the simulated trajectory, the parameter set was also required to retain non-negativity of each state 1 ( , ) , and first-order parameters (Equation 5). Each of these functions is a simple penalty function which finds the most negative value that violates the constraints, and calculates a weighted value that is added to f , to further constraint the parameter set, p v (gives a number between 0 and 1, which can be further weighted by α and β ) in this optimization procedure. Moreover, the algorithm was run for 300 function evaluations (which on average was roughly 100 iterations past convergence) to ensure that each model parameter set had reached convergence. The primary reason for using this global optimization method was to avoid the potential systematic error of landing in local minima repeatedly, thereby severely biasing the results.

Multivariate sensitivity analysis
Once the parameter set defining the best fit to the experimental data for a particular network model ('fitted set') was determined, we performed multivariate sensitivity analysis based on a percent change of the all the parameters in the fitted set varied simultaneously. Since we were interested in the local control properties for each parameter, only small, local perturbation of 5% of the nominal parameters value (obtained from fitting) were tested. The purpose was not to find bifurcation boundaries or other features of the global system behavior, as we were not interested in cases where the simulated dynamics were dissimilar to the experimentally observed behavior. The only criterion was that the perturbation be sufficient to determine sensitivity of a given model to variation of different parameters. Thus, for each model we compared the results of 200 distinct implementations, each time sampling the parameter values from a uniform distribution contained within +/-5 % of the nominal value found in the fitted set. In the event that a parameter set broke any of the aforementioned constraints, the parameter set was discarded, and another simulated in its place. For each parameter generation, the resultant trajectories describing the proportions of cells in each state were compared to the corresponding experimental data, and the sum of squared residuals (SSR) value was computed for all four curves representing the evolution of all four states (Supplementary Supplementary Figure 16). The results were assembled for all 298 topologies to be used in subsequent model ranking and parameter control analysis (Section 3, below).
In addition, for each of the 200 simulations, and 200 corresponding SSR values, a correlation coefficient was calculated for the values of each parameter and the corresponding SSR metric. To ensure that the statistics didn't reflect unbounded trajectories, each trajectory was checked for non-negativity, both in the state dynamics x, and in the state transitions, h. In the event that a simulated trajectory or its transition was negative, that parameter set was discarded, and another was generated. At each iteration, the parameters and statistics are saved into a list of object struct in their respective fields (seen below).
Nodes(i).field At each iteration, the parameters and statistics are saved into a list of object struct in their respective fields (seen below).
Nodes(i).field  plotted as a histogram (see e.g., Supplementary Fig. 18). The peak of the histogram corresponds to the most common goodness of fit yielded by a particular model topology. The variance of the SSR value distribution given by the histogram reflects the sensitivity of the goodness of fit to parameter variation. The models with relatively low peak histogram values are considered to give better fits to the experimental data. The models with relatively low histogram variance values are considered to display robustness of the fit to parameter variation. Since one knows for all model topologies whether or not a specific interaction between states is present or not, one can ask the question: do certain parameters occur more frequently in models with better fitness and/or high robustness to parametric variation? If a quantitative metric giving an answer to this question is developed, one can rank-order the parameters according to the values of this metric and postulate the relative likelihood that a certain cell-cell interaction is present in the developmental system under investigation.
We thus used the properties of the histograms of the types shown in Supplementary Fig. 18 to define the metrics.
The goodness of fit for a particular model was defined as the position of the peak of the histogram. Then for each parameter, we analyzed the goodness of fit of the models, in which the interaction described by this parameter is present.
Since we wanted to give higher weights in the rank ordering to parameters occurring in models providing better fits, we have calculated for each parameter the following metric: where fM is the mean of the distribution represented by the histogram of the SSR values for M-th model, for the j-th parameter of interest.
Another way of rank-ordering model parameters is to segregate the parameters, which appear to endow resistance to parametric uncertainty or robustness. The underlying assumption is that if a model should not only fit the data well, but also do so robustly, so that the fit is not easily altered if the parameters defining the model are slightly altered. Thus we again considered for each parameter j the models, in which this parameter is present. We further calculated the robustness metric Rj, which is the summation of the inverse variance values of the SSR value histograms corresponding to each of these models M.
Therefore, the robustness metric value is greatest for the parameters participating in the models with fits generally resistant to parametric uncertainty. We have rank-ordered the parameters according to the metrics F and S and plotted the results in Supplementary Fig. 19. We have represented the relative input of models, in which the values of the corresponding parameters are either positive or negative, into the sums defining the metrics, by using different colors. Thus one can see whether the signs of the corresponding interactions are either positive or negative. One can see that 1) there is a good correspondence between two ways to rank-order the parameters and 2) the robustness metric has greater variability allowing for a better discrimination of the relatively important parameters. Thus we have chosen to primarily focus on the robustness metric in our further analysis.
A note on the analysis: Although it might appear that ranking models based on SSR statistics alone ignores the fact that different models have different numbers of parameters and therefore this analysis should be normalized by the number of parameters, as done in e.g., Akaike or Bayesian Information Criteria, this is in fact not true. We stress that our analysis does not rank models, but rather, it ranks interactions (the network links). For each such interaction, we calculate the metrics Rj and Fj that rely on the performance of the entire set of all models in which it is present. Because any given interaction participates in an equal number of models of any given complexity (that is, models that contain one, two or three interactions), the process of ranking interactions described below is properly normalized. Specifically, each of the 12 possible feedforward and feedbaclk interactions participates in 67 models, of which 1 contains that interaction alone, 11 contain that interaction together with one other, and 55 (=11x10/2) contain this interaction together with two other. Thus, although the analysis of specific models is uncorrected for the number of parameters, the analysis of the entire set of models containing a specific interaction is unbiased: each relies on the same number of models of any given complexity.  Fig. 20. From Fig. 20 it also becomes more evident that ff23, ff34 and fb43 are present more frequently in more robust models, whereas fb32 is present mostly in the least robust models. The rest of the parameters show mixed distributions and often mixed sign in Supplementary Fig. 20, making any judgment about their putative roles in the underlying erythropoiesis network difficult. These preliminary results were tested by determining their sensitivity to alterations of either the experimental dataset or the methods of fitting and analysis (Section 4). The final discussion of results was made following this further analysis.
Since each fitting procedure also results in specific values of the parameters in various models, we have plotted the values of the parameters to determine if any characteristic values emerge ( Supplementary Fig. 21). The interesting features of this dataset include relative magnitude of transition rates: p23 > p12 > p34. Additionally, it shows relatively high magnitudes of fb43, ff23 and fb32. As noted above, fb32 occurs mostly in models of low rank and thus probably does not occur in the true erythropoiesis regulatory network. The results show that the metric R profiles display little sensitivity to this variation. Therefore it is likely that although some properties of the distribution of R may alter, the relative rankings of ff23, ff34, fb43 (all high) and fb32 (low) are not sensitive to this data set perturbation. now ranked as less significant. Importantly, however, ff23 is again ranked as the most significant with the robustness metric 0.3, further suggesting the importance of this interaction. The control diagram also reveals that ff34 and fb43 lose their control properties in highly ranked models, but ff23 retains a significant amount of control.

Removing E14.5, 2.5 days to full differentiation
In order to complement the analysis in section 4.2, where E15.0 and E15.5 were kept constant with respect to time, we also assumed that the entire interval, the differentiation terminating with Ter119 HIGH , CD71 LOW occurs over 2.5 days, as opposed to 3.0 days. The results show that although now ff23 is still ranked topmost in terms of its significance, its robustness metric value, 0.3, is more closely matched by many other possible interaction mechanisms, so that the overall spread of robustness metric values (except fb32 and ff24) is very narrow, making the ranking results less reliable.
Nevertheless, both ff23 and ff34 are ranked at the top again.

Testing Computational Methods
We also studied the sensitivity of the results to exact computational methods of determining the predominant model topologies. In alternative analyses, we studied suitability of relying on particular criteria for model selection. In the algorithm we used, data on fractional cell numbers in all four states were used in fitting a model to the experimental data and estimation of residuals of the fitted model and 200 corresponding model perturbations. However, the precise gating of the first and last stages is subject to more uncertainty compared to the stages 2 and 3, based on the magnitude of ratios of cells in these stages in each of the experimental data points. The first stage and the last stage primarily occur in the first and last time point measurements respectively, and as the proportions of cells in these states tend to be low, there is a greater possibility of error, e.g., due to minor changes in gating. Therefore, as an alternative method, we tested the role of increasing reliance in treating regression in some cases, and sensitivity in others, on the middle two stages, which appear predominantly throughout the course of embryonic development. Due to propagation of the erythropoietic wave, comparatively less data exist on the dynamics of cell numbers in states 1 and 4. We thus tested if the results would be affected if the analysis relied primarily on states 2 and 3. As the system is conserved, every parameter is represented even in this reduced form, as the differentiation equations for stages 2 and 3 contain all transitions, h12, h23, and h34 (Equations 1 and 2). Therefore, parameter optimization still reflects the first and last stage, but is weighted much less in their favor. As can be seen in Supplementary Fig. 26, the two predominant mechanisms become even more dominant, and are highly favored as the regulatory connections in this system.
Supplementary Figure 26. Ranking of interactions according to the R metric and the corresponding control matrix with states 2 State Regression, 2 State sensitivity, 50% parameter variation We initially chose to test only local perturbations of parameter values (5%). To test if the previous results hold at very large perturbation we ran a simulation where we used a 50% parameter variation. Supplementary Fig. 27 shows that even in this case, ff23, ff34 and fb43 are top ranked.
Supplementary Figure 27. Same as in Supp. Fig. 26, but 50% variation of parameters in multi-variate sensitivity analysis 2 State Regression, 2 State Sensitivity, similarity to the fitted data In the preceding analysis, in parameter variation studies, we compared each generated time course to the experimental data given in Table 1. It is of interest to see if the results would change, if the comparison was instead performed with respect to the curves resulting from the fitting procedure for each model. The results are shown in Supplementary Fig. 28.
The ff23 term remains ranked at the top, but the rest of the results are significantly less differentiated. We note that this analysis does not reflect the experimental data to the same degree as the analysis above, since we measure the similarity during the sensitivity analysis between the solutions resulting from the perturbed parameter set and the original fit, so the similarity to the experimental data is confined to the fitted set. The parameters that appear robust in this analysis may not actually endow the system with a good fit. were used in fitting and parameter variation. As can be seen in Supplementary Fig. 28, ff23 is still top ranked and fb43 is the second top ranked parameter, although the value of S for fb43 is similar to those of most other parameters.
Supplementary Figure 29. Same as in Suppl. Fig. 28, but 4 states are used in fitting and multi-variate sensitivity analysis 1.6. Relative rankings of models with ff23 and ff34 interactions As the main predictions of the model were that ff23 and ff34 would be likely to be found in any model describing fetal erythropoiesis, we wanted to see how the models containing these two interactiosn would rank compared to other models. However, we need to stress that the compatibility of any two interactions (including ff23 and ff34) is not in and of itself a function of the algorithm's performance. Each interaction participates in an identical number of models for each given model complexity (i.e. models with one, two or three interactions). Of the 298 model topologies we tested, the ff23 and ff34 interactions co-exist in one model where there are only two interactions, and in 10 models in which there are 3 interactions. As can be seen in the Supplementary Figure 30 we find that these models are indeed ranked among the best models in terms of both their fitness and robustness. We also confirmed through Wicoxon non-parameteric test that both the fitness and robustness of models with ff23,ff34 interactions are significantly better (p<0.001) than those in the rest of the model population.
Supplementary Figure 30. Ranking of models by the mean and variance of the SSR (left column, SSR is equivalent to 'Residual Sum'), and box plots contrasting the distributions of SSR means and variances for the entire model population and the subset of models containing both ff23 and ff34 (right column). The blue circles indicate models containing both ff23 and ff34, whereas the red circles indicate the model containing ff23 only. The rankings of the two models with just ff23 or just ff23-ff34 are indicated. The box indicates the inner quartile range. The median is the horizontal red line within the box. The whiskers extend to cover the upper and lower quartiles up to a distance of 1.5 times the inner quartile range. Red crosses indicate extreme values (outliers beyond the whisker length).

Interpretation of prominence of inter-state interactions
It is clear from the Supplementary Figs. 19-29 that the analysis of both the original and modified (control) data sets consistently yields similar results in terms of the parameters deemed significant for maintaining robust fit of a model to experimental data. We hypothesized that cell-cell interactions represented by these parameters might exist in erythopoiesis.
In particular, three interactions stand out as ranked at the top with the highest significance across the control simulations.
These are the positive feedforward interactions ff23 and ff34 and the negative feedback interactions fb43. The first two interactions represent (see Supplementary Table 3) increased rates of transition out of the states 2 and 3 in proportion to the squares of the numbers of cells in 2 and 3 respectively. Thus these mechanisms are autoregulatory in the sense that transition out of a particular state depends only on the population in that stage; both the source and target for regulation are the same. Since both ff23 and ff34 are predicted to be positive, their effect is to decrease the number of cells in states 2 and 3, in proportion to the relative numbers of cells in each of the states, indicating a clear possibility that either an increase of the number of cells in these states increases the rates of differentiation out of these states or the rates of apoptosis in these states. In the latter case, some signal must exist that would allow the cells in either state 2 or state 3 to act in autocrine or paracrine fashion on the cells in the same state leading to apoptotic death. We hypothesized that that this autocrine and paracrine pro-apoptotic mechanism indeed is responsible for the observed behavior, focusing on Fas and FasL interactions as potential candidates. As described in the text, the experimental results provide support for this hypothesis.
The negative feedback interaction fb43 corresponds to an effective decrease of rate of transition between states 3 and 4, in proportion to the product of the relative numbers of cells in these states. This feedback can be interpreted as either upregulation of the number of cells in the state 3 due to a decrease rate of either differentiation or apoptosis from this state (both as a function of the number of cells in the state 4), or as an increased cell death in the state 4, in proportion to the number of cells in the state 3. The latter interpretation is a consequence of operating with normalized cell numbers and thus conservation assumption made for the model, and of the absence of a transition out of the state 4. Thus, if the cells in the state 4 undergo apoptosis, the effect will be in apparently less efficient transition from 3 to 4, which is the source of populating the state 4. The data provided in the main text agree with this latter interpretation, although the former ones cannot be excluded.
Interestingly, the parameters responsible for the interactions ff23, ff34 and fb43 are either predicted to be strongly positive (the first two) or negative (the last) unlike many other parameters, whose sign may vary with the actual model ( Supplementary Fig. 21). In addition, ff23, ff34 and fb43 are predicted to exercise a substantial degree of control of the fit as predicted by the fact that they occur with high frequency in highly robust models and that the values of the control metric for the corresponding parameters in highly robust models are relatively high ( Supplementary Fig. 20 and the right panels in Supplementary Figs. 23-29). These observations further support potential importance of these interactions and argues for their presence in the real erythropoietic system. Below we provide the proposed scheme of interactions in the real erythropoiesis developmental regulatory network. In addition, to including interactions ff23, ff34 and fb43, we have assigned different thickness to unmodulated transition rates represented by the parameters p12, p23, p34, as determined from Supplementary Fig. 21.