Cell fate clusters in ICM organoids arise from cell fate heredity and division: a modelling approach

During the mammalian preimplantation phase, cells undergo two subsequent cell fate decisions. During the first decision, the trophectoderm and the inner cell mass are formed. Subsequently, the inner cell mass segregates into the epiblast and the primitive endoderm. Inner cell mass organoids represent an experimental model system, mimicking the second cell fate decision. It has been shown that cells of the same fate tend to cluster stronger than expected for random cell fate decisions. Three major processes are hypothesised to contribute to the cell fate arrangements: (1) chemical signalling; (2) cell sorting; and (3) cell proliferation. In order to quantify the influence of cell proliferation on the observed cell lineage type clustering, we developed an agent-based model accounting for mechanical cell–cell interaction, i.e. adhesion and repulsion, cell division, stochastic cell fate decision and cell fate heredity. The model supports the hypothesis that initial cell fate acquisition is a stochastically driven process, taking place in the early development of inner cell mass organoids. Further, we show that the observed neighbourhood structures can emerge solely due to cell fate heredity during cell division.

. Flowchart of the implemented model.

2/11
Parameter estimation for different hypotheses The model can be used to test several assumptions addressing cell fate heredity. In total, four different hypotheses are tested. Each hypothesis is considering cell fate switches during cell division according to Fig. 1d. During cell division, the cell fate is passed on to the daughter cells. A cell fate switch is possible with a given rate. The model considers the following cell fate switches: N -Gremain N -G -(ζ 1 ) or become N + G + (α 3 ). N + G + remain N + G + (ζ 2 ) or become N + G -(α 1 ) or N -G + (α 2 ). N + Gand N -G + remain N + G -(ζ 3 ) and N -G + (ζ 4 ) or switch to the opposite cell fate (α 4 ) and (α 5 ), respectively. These cell fate transitions form a system of linear ordinary differential equations, which can be written as with the analytical solution: with v and λ the eigenvectors and eigenvalues of the coefficient-matrix A, respectively. The unknown prefactor c values can be determined by inserting the known cell counts of the different cell types at t 1 . The rates for the different cell fate switches vary for the different hypotheses (see Tab. A1).

Hypothesis 2, 3 and 4
The neighbourhood statistics for H 2 agreed reasonably well with the experimental data (see Fig. A3). The neighbourhood pattern of cells with the expression types N -Gand N + Gat t 1 were similar to the neighbourhood structures obtained under assumption H 1 , including the misfit involving N -Gcells. In addition, simulated N + G + cells were significantly less often neighbours of other N + G + cells than the experimental data suggests (see Tab. A4). The latter speaks against a cell fate change of N + G + cells into N -G + cells. If performed under the assumption of the H 2 model, the evaluated proporitons for t 2 showed a better agreement with the experimental data. In particular, the predicted neighbourhood statistics of N + G + cells differ statistically less substantially from the in vitro measured proportions compared to the simulated cell neighbourhood statistics under the assumption H 1 (see. Fig. A2 and Tab. A2).
If we consider cell fate switches from N + Gto N -G + and vice versa, which according to the literature are considered as unlikely, a third and a fourth hypothesis can be formulated. In both hypotheses, the strong increase in the amount of N -G + is explained by cell fate switches from N + Gto N -G + . While H 3 permits only the cell fate switches from N + Gto N -G + , H 4 is considering a small flux between both cell fates. The parameter values for both hypotheses are shown in Tab. A1. As expected, the simulated proportions for H 3 and H 4 agree very well with the experimental data at t 1 and t 2 (see. Fig. A2 and Tab. A2).
Comparing H 3 and H 4 to H 1 and H 2 , we observe that their ψ values are higher, thus the quality of the fit of the neighbourhood statistics is lower (see.

8/11
p-values Table A2. p-values for the statistical comparison of the simulated expression type composition of ICM spheroids to the experimental data at t 1 and t 2 . Simulations were performed under the assumption H 1 , the assumption H 2 , the assumption H 3 and the assumption H 4 . For statistical comparison a Wilcoxon-Mann-Whitney test with Bonferroni correction is performed. Statistically significant differences between simulated and experimental data are indicated in red (p < 0.05).  Table A3. p-values for the statistical comparison of the simulated expression type composition of neighbouring cells of ICM spheroids to the experimental data at t 1 and t 2 . Simulations were performed under the assumption H 1 . For statistical comparison a Wilcoxon-Mann-Whitney test with Bonferroni correction is performed. Statistically significant differences between simulated and experimental data are indicated in red (p < 0.05). Not statistically significant differences between experimental data and simulated expression type compositions of neighbouring cells for H 1 (excluding N -Gcells) with t 0 = 200 cells are indicated in blue.
neighbouring cell t 1 t 2 expression count expression type expression type type at t 0 N -G -   Table A6. p-values for the statistical comparison of the simulated expression type composition of neighbouring cells of ICM spheroids to the experimental data at t 1 and t 2 . Simulations were performed under the assumption H 4 . For statistical comparison a Wilcoxon-Mann-Whitney test with Bonferroni correction is performed. Statistically significant differences between simulated and experimental data are indicated in red (p < 0.05).