The model is initiated by a population of cells, the properties of which are evolved by specific pre-assigned rules. The cells are primarily categorized into two signaling regimes: Ω and A; Ω can be considered as an active regime supporting cellular proliferation and differentiation, while A is a more dormant regime where cells are quiescent and prone to dedifferentiation. The cells can transfer in between these two regimes, an event decided upon primarily by a cell-specific parameter, termed as ‘a’ value. This parameter ‘a’ is randomly assigned to each cell at the beginning of simulation, and is updated at each time step. The probability of transfer to/from a regime is dependent on this parameter along with the number of cells in the destination regime. This ‘a’ value is unaltered in the A regime while it progressively reduces in the Ω regime, and when it falls below a threshold the cells lose their ability to transfer to the A regime.

Each cell is randomly assigned a maximum life span, exceeding which it will die. While the cells age in the Ω regime, they neither proliferate nor age in the A regime. Proliferation is allowed in the Ω regime for an amount of time which is cell dependent, after which the cell enters a senescent stage and will not proliferate. Furthermore, an individual cell is allowed to proliferate only after it loses the capability to pass into A regime having crossed the ‘a’ threshold value.

Cellular differentiation is governed by the ‘lineage propensity’ parameter, representing a cell's likelihood to differentiate into a particular lineage. Only the Ω regime allows an increase in lineage propensity. While updating the propensity of differentiation to a particular lineage, all the possible lineages are competing and any can be updated, with the one with higher propensity having a higher probability of being selected. Once a cell's propensity exceeds a threshold level, the cell is considered committed to that particular lineage and will retain its differentiated phenotype. If this threshold has not been exceeded and if the cell is chosen to be transferred to the A regime, the propensity values will converge to an average value. The model is therefore able to track specific germ layer populations, and through this the percent of the population positive for Sox17 (visceral and definitive endoderm marker) and CXCR4 (definitive endoderm and mesendoderm marker) [17]–[19] can be extracted.

The current work is focused on the mechanistic investigation of the dynamics of hESC induction into endoderm. Using the platform of the stochastic population based model we investigated several alternate mechanisms and analyzed them for agreement with experimental data. Three characteristics of the differentiation process were chosen to be investigated: the presence/absence of an intermediate germ layer, mesendoderm, which subsequently gives rise to mesoderm and endoderm [20]; the presence/absence of CXCR4, in mesoderm ; and whether proliferation of a specific differentiated cell phenotype is favored over others. Combination of aforementioned attributes results in 12 alternate mechanisms, each of which were incorporated into the model and analyzed for agreement with experimental data. It is expected that the most likely mechanism will best describe the experimental dynamics of the stem cell system.

The incorporation of mesendoderm involved a two stage differentiation scheme. In the first stage, hESC are able to differentiate into either mesendoderm or visceral endoderm. Once cells are committed to the mesendoderm lineage, several of their attributes are re-initialized, such as their ‘a’ value and lineage propensities. The mesendodermal cells can then further differentiate into endoderm or mesoderm. Differences in the proliferation potential of different phenotypes were incorporated by considering 3 scenarios: proliferation of hESC and mesendoderm; proliferation of hESC, mesendoderm and endoderm; and proliferation of all phenotypes.

As with any stochastic model, the number of model runs necessary to obtain a converged solution needs to be determined. An additional parameter of the current model is the initial cell population, which affects the solution over a certain range. A two-dimensional convergence study was undertaken, wherein the effects of both stochastic run number and initial cell population on model output were determined. The convergence test allows determination of the minimum number of stochastic runs and initial cell population beyond which the model output does not significantly change. All results reported here are using the converged parameter values.

Sensitivity analysis was performed to determine the relative importance of parameters in affecting the outputs of cellular growth, death, and lineage commitment. Because this model is probabilistic in nature, traditional ways of determining local sensitivity, e.g. partial derivative of output with respect to an input, cannot be employed. A stochastic analysis was therefore chosen, using differences in output histograms of nominal and perturbed parameters, S, to estimate parameter sensitivity [21] as:(1)where x_j and y_j are the individual elements in the nominal and perturbed histograms, respectively, I_i represents the range of each bin i in the nominal histogram, X is a counting variable which takes on a value of 1 if x_j/y_j is within the interval I_i, |x| and |y| are cardinalities of data sets x (nominal) and y (perturbed), respectively, and k is the number of bins, which was determined by calculating the appropriate bin size of the nominal output histogram by the Freedman-Daconis rule [22]:(2)Where IQR(P) is the interquartile range of a sample population P and n is the number of observations in P. The number of bins was different for each output, and ranged from 37 to 51.

For sensitivity analysis, each parameter was perturbed by 10% while keeping the rest at the nominal value. For each bin in the nominal output histogram, the difference between the percentage of total nominal histogram elements residing in that bin and the percentage of the total perturbed histogram elements residing in that bin is calculated. The sum of the absolute value of this difference over all of the bins is the histogram distance.

Having determined the most sensitive model parameters, the next step is to determine an appropriate value of the parameters which best estimates the experimental data. We realize that a single parameter combination may not be adequate in describing the experimental data; instead, there exists a parameter hyper-space adequately satisfying the data. Hence an ensemble parameter estimation was performed by randomly generating initial guesses from the hyper-space of the sensitive parameters. The model was simulated with 10000 random parameter samples; for each of these simulations the least square estimate is determined between experimental data and model output. These simulations were run for each mechanism and condition under investigation, and parameter ensembles were generated by considering only those parameter sets which meet certain error constraints. Following the above detailed methodology we evaluated the model predictions obtained from the different mechanisms and compared them with experimental data to determine the most plausible mechanism.

In the next step the developed stochastic model was used to test the proposed mechanisms for agreement with experimental data. The mathematical model involves multiple parameters which require detailed analysis before the model can be used for prediction. The parameters can be grouped into two categories: (a) simulation parameters which affect the convergence behavior of the simulation and (b) model parameters which affect specific model output for the converged simulation. Two parameters were identified to be simulation parameters: initial cell population and number of stochastic runs. These parameter values were optimized by performing a two-dimensional convergence study, as illustrated in Figure 2 for the CXCR4-positive population output. Overall it is observed that the initial cell population has a more dominant effect on convergence, while the effect of stochastic runs was rather weak beyond 2000 runs. Following this analysis an initial cell population of 9000 and total stochastic runs of 4000 was used for subsequent simulations.

A detailed sensitivity analysis was performed for the model parameters in order to determine the relative importance of the parameters in determining model output. As detailed in Equation 1, the measure of histogram distance is used to represent the parameter sensitivity associated with a specific model output. Figure 3A illustrates the model parameter sensitivity to the output of cellular growth, as concluded from the shift in histogram distance (inset). A clear jump in the sensitivity is observed, with a large difference between parameters with low and high sensitivity.

While Figure 3A represents the parameter sensitivity to model output of cellular growth, similar analysis was also performed for all of the model outputs (data not shown). Overall it was observed that even though the magnitude of sensitivity differs between outputs, the highly sensitive parameters were mostly conserved across outputs. Furthermore, in the present study we are investigating multiple competing mechanisms, which require modification of the model formulation. Since the effect of such modifications on the parameter sensitivity is not intuitively obvious, similar analysis was repeated for each of the 12 proposed mechanisms. Figure 3B summarizes the number of sensitive parameters for each of the mechanisms. From the analysis, eight classes of parameters were consistently observed to have the highest sensitivity:

Having determined the sensitive parameters for each of the mechanisms, the next step is to determine the optimum parameter values which will result in best agreement with experimental data. In literature, biological samples are typically defined as being ‘sloppy’ [23] with a broad ensemble of parameters satisfying the error constraints. Accordingly we also target identification of representative ensemble of parameters. The model is formulated to capture the dynamics of cellular growth, death and differentiation, the output of differentiation being of most interest. Hence the model parameters were optimized with respect to differentiation dynamics, while growth kinetics and the dynamics of cell death were used for verification. A projection of the simulated error onto a 2-dimensional parameter space (for the mechanism which incorporates mesendoderm and promotes proliferation of both uncommitted and endodermal cells without CXCR4 being expressed in mesoderm (‘Mechanism B’)) is shown in Figure 4A. Although it was initially thought that a trend might be observed between the error and values of the parameter ensemble, Figure 4A shows that there is a lack of any correlation between multiple parameters and associated errors (shown for parameter ‘d’; further analysis of all parameter combinations yielded same results, data not shown). Figure 4B illustrates the minimum ensemble error for each of the proposed mechanisms simulated under the two endoderm induction conditions, the error being evaluated according to least square formulation.

As shown in Figure 4B, the absence of mesendoderm gives rise to large errors, in some cases an order of magnitude higher than their counterpart models which include mesendoderm. If one considers the Activin A only condition, the most accurate mechanisms include those which incorporate mesendoderm and promote proliferation of both uncommitted and endoderm germ layer both with (‘Mechanism A’) and without (‘Mechanism B’) CXCR4 in mesoderm. Since Figure 4B illustrates the accuracy of the model in predicting differentiation dynamics only, the performance of the 2 prospective mechanisms were further verified with the help of growth kinetics and the dynamics of cell death. Figure 5 illustrates the ensemble simulation of all the model outputs and its comparison with experimental data. While both the mechanisms had excellent performance in predicting Sox17 and CXCR4 dynamics (Figure 5 C,D,G,H) they differed significantly in predicting growth kinetics and cell death dynamics (Figure 5 A,B,E,F). Figure 5 clearly illustrates that Mechanism B performs better in describing both growth kinetics and cell death dynamics compared to Mechanism A. Hence the former was chosen to be the most likely mechanism for Condition A.

Condition B (Activin A supplemented with FGF2 and BMP4) proved more difficult to describe via the investigated mechanisms, mainly because CXCR4 dynamics exhibits a faster and more prominent drop as compared to Activin A only condition. As shown in Figure 4B, the two mechanisms which give lowest error for Condition B are the ones which incorporate mesendoderm, have CXCR4 present in mesoderm and promote proliferation of all phenotypes (‘Mechanism C’) and the previously described Mechanism B. The simulated dynamics of these two mechanisms with experimental data of Condition B are shown in Figure 6. As with Condition A, although the incorporation of CXCR4 in mesoderm gives a small error, the simulated profiles of growth kinetics and cell death are not in agreement with experimental data. The next best mechanism, Mechanism B, exhibits both a low error and good results with all outputs, so it was again chosen as the most likely mechanism.

It is therefore reasonable to conclude that during endoderm induction with the conditions described above, undifferentiated stem cells first differentiate into a mesendoderm germ layer with subsequent differentiation to endoderm and mesoderm, the latter not expressing CXCR4. Furthermore, the induction condition seems to promote proliferation of both pluripotent and endoderm-like cells. The optimized parameters of this mechanism are shown in Table S2, with definitions of parameters in Table S1.

The power of mathematical models lies in their predictive capacity. The predictive capacity of our proposed model was thus tested by simulating the population dynamics of cell types for which no a priori data was used in constructing the model. The chosen populations were that of undifferentiated cells and mesendoderm cells. The simulated profile of the undifferentiated cells (Figure 7A and B) shows an exponential decay to a final value of 10% of the cellular population. This final value was reached in approximately 3 days. The mesendoderm cell population was predicted to display more interesting dynamics, with a transient increase in cell population over the first day, followed by a decreasing trend over the next few days (Figure 7D and E).

In the next step the validity of such prediction was verified by conducting further experiments to analyze the dynamics of undifferentiated cells by Oct4 gene expression and that of mesendoderm cells by Brachyury expression. While the comparison of population dynamics with mRNA levels is not exact, under the assumption of efficient translation they become comparable. Figure 7 illustrates the comparison of experimental data with model prediction, which are found to have excellent agreement given that the model was generated with no information of these specific cellular dynamics. Oct4 levels exhibit a decay to a final value of around 20% of the maximum, at a time which correlates with simulated predictions (3 days). Brachyury levels showed a similar bimodal trend as was predicted by the model. It reached a maximum around 24 hrs, following which it gradually decayed over time.