The raw data is available at GEO GSE14782.

The processed adrenal.human.seurat.scrublet.rds dataset used in this study can be found at http://pklab.med.harvard.edu/artem/adrenal/data/Seurat/.

This represents single-cell transcriptomic analysis using 10x chromium on isolated individual cells from dissected adrenal glands with surrounding tissue from human embryos at 6, 8, 9, 11, 12, and 14. The initial matrix displays the expression of 25787 genes across 72571 cells. Only cells harboring a cell doublet probability (Scrublet score) <0.2 were kept.

The cell types were annotated by the authors based upon gene expression pattern.

To elucidate the sympathoadrenal cell fate transitions in humans at a higher resolution, we focused on the 3,901 cells of SCP, chromaffin, and sympathetic fates, omitting annotated cell cycle genes. As advocated by the authors, we removed a small potentially contaminating population of cells expressing the STAR gene. This resulted in a final matrix displaying the expression level of 20807 genes in 3759 cells at six developmental time points. We verified that the UMAP generated from this dataset was in perfect agreement with Figure 2a in the original paper (see Figure 1a). UMAP was computed using the Uwot package. ²² All the UMAPs were embedded within the experimental UMAP using a combination of the ret_model and umap_transform functions of the uwot package.

The stemness score was computed per cell as signature_SCP/signature_SCP + signature_sympatoblast + signature_chromaffin, by summing up the gene expression level for the following signature genes: -

SCPgenes: PLP1, FOXD3, FABP7, S100B, NGFR, ERBB3, MBP, MPZ, COL2A1, POSTN, MOXD1, and GAS7

First, we cloned the original CARDAMOM git repository at https://github.com/ogandril/cardamom. All computations were performed at the IN2P3 computing center ( https://cc.in2p3.fr). The scripts used for launching CARDAMOM on a distant server are available at https://github.com/ogandril/Launching_Cardamom. All interaction values were scaled by a factor of 10 and thresholded at a value of four to maintain only the most relevant interactions.

The resulting inter.csv interaction matrix was opened with R ²³ and plotted using the qgraph package. ²⁴ A list of transcription factors was obtained from https://github.com/saezlab/CollecTRI. ²⁵ Violin plots were drawn using the ggbeeswarm package.

The simulation was performed using the CARDAMOM build-in version of HARISSA, a software for performing GRN simulations based on an underlying stochastic dynamical model driven by the transcriptional bursting phenomenon. ²⁶

The dynamics of each gene is given by the following ordinary differential equations: M ` i = s 0 , i E i − d 0 , i M i P ` i = s 1 , i M i − d 1 , i P i where M _i is the number of mRNA molecules for gene i, P _i is the number of proteins for gene I, E _i is the state of the promoter, s _0,i is the mRNA synthesis rate, d _0,i is the mRNA degradation rate, s _1,i is the protein synthesis rate, and d _1,i is the protein degradation rate.

E _i can randomly switch from 0 to 1 (activation) with a protein-dependent rate k _on,i (P) and from 1 to 0 (inactivation) with a constant rate k _off,i.

We define the protein-dependent rate function k _on,i as in Ref. 26. k on , i ( P ) = k min , i + k max , i exp ( β i + ∑ j θ j , i P j ) 1 + exp ( β i + ∑ j θ j , i P j ) such that the promoter activation rate is between k _min,i and k _max,i. The parameter β _i represents the basal activity of gene i, and the parameter θ _j;i encodes the value of the j → i interaction (see below for its dertemination).

Because we simulated a mechanistic model, one needs to provide HARISSA with half-lives for both mRNAs ( d _0,i ) and proteins ( d _1,i ). Experimentally determined values for mRNAs and proteins were extracted from ²⁷ and, ²⁸ respectively. Genes that were not found in these datasets were attributed to the mean value of experimentally determined half-lives. Degradation rates were obtained using the following formula: degradation rate = ln(2)/half-life. The half-lives of proteins were limited by a parameter called the cell cycle, which was set to 20h.

Color code for a rapid estimation of the fit quality was as follows: given all Kantorovich distances for all genes at all times, we considered that 40 % of the values should be considered as correct (after a careful examination of the marginals) and therefore colored in green. The remaining 60% was split into two, the first half colored in orange, and the last half colored red.

In CARDAMOM,we compute at each time point t _i with i ≥ 1 the value of each interaction, using θ i ≔ argmin { θ ∑ c = 1 N i ∥ k on θ ( P c ) − k P c ∥ 2 } + λ ∥ θ − θ i − 1 ∥ 2 where Ni is the number of cells at time t _i, P _c the amount of protein in cell c, k P c the frequency mode and θ 0 = 0 .

So we do obtain at each time t _i a variation in the network ∆ i ≔ | θ i − θ i − 1 | , which quantifies how much each interaction estimate has been impacted by the cells measures at that time compared to all previous times. We can therefore, for any interaction between genes k and l, computes the times for which the largest impact was observed using t k , l ≔ argmax i = 1 , … , T ∆ i k , l .

T being the latest time point available.

The dynamical interactions Θ i k , l can therefore be defined as: Θ i k , l = θ T k , l 1 t k , l = i .

1 being the indicator function.

Our first task was to organize the cells according to their differentiation status, as our GRN inference is anchored in the time-dependent evolution of the differentiation process. Cells were collected from human embryos at 6, 8, 9, 11, 12, and 14 weeks post-conception. First, we assessed whether this cell collection time could be a relevant framework for generating time-reconstructed data. As shown in Figure 1b, this was not the case, since all cell types were present at all collection times. This is in line with the fact that in the human adrenal gland, the transition from Schwann cell precursor ( SCP) to sympathoadrenal states persists for several weeks, ¹⁰ which leads to the expectation that cells at all stages of differentiation should be present at all collection times.

Therefore, we were faced with a different situation from the one we used to benchmark CARDAMOM, where the cell sampling time was in accordance with the actual process time. ¹⁸

Consequently, we organized the cells according to their stemness scores ( Figure 2a).

For this, we assumed a differentiation sequence starting with the SCP, giving rise to chromaffin cells and sympathoblasts through a bifurcation. Such a scheme is in line with the conclusions of the original paper from which we obtained dataset ¹⁰ and is shared by most ^{6, 12, 13, 29} but not all ^{4, 30} authors working on sympathoadrenal differentiation in humans. This led to the definition of the stemness score (see Materials and Methods) as the ratio of SCP-signature gene expression, evolving from 1 (only the SCP signature genes are expressed) to zero (none of the SCP signature genes are expressed). This allowed us to define the “process time” ( Figure 2b) using an overlapping sliding window of 300 cells. This was defined only on the first 1300 cells since the stemness score was almost down to zero at the end of that sequence, and because we wanted to capture in detail the abrupt decrease in the stemness score. The use of an overlapping window is justified by the fact that differentiation should be considered a discontinuous nonlinear process (see e.g., Ref. 31), and therefore should allow cells to go back and forth in the gene expression space. This resulted in a 20807 genes x 1800 cells matrix.

We then explored the most relevant gene set to be used for the GRN inference. Because our inference and simulation scheme CARDAMOM is fundamentally rooted in exploiting the power of distributions, ^{18, 32} we first computed a Kantorovitch distance ³³ for all genes for all delta-times. Similarly, since most delta-entropic genes were shown to represent the most relevant genes for the differentiation process, ³⁴ we also computed delta-entropy for all genes for all delta-times.

The combination of these two lists ( Figure 3) provided us with a list of 97 genes ( Table 1) that were used for inferring the GRN.

This left us with a matrix displaying the expression levels of 97 genes in 1800 cells distributed over six time points. We verified that the information contained within this matrix was sufficient to correctly capture the differentiation process of interest, as assessed by UMAP representation ( Figure 1c). It is important to note that these genes were chosen purely based on their expression patterns, irrespective of their putative biological functions (see Discussion).

We then applied the CARDAMOM algorithm to infer the GRN structure ( Figure 4a), resulting in a relatively densely connected network.

First, we explored the network connectivity. We observed a striking difference between the distribution of the numbers of outgoing and incoming nodes ( Figure 4b). The incoming nodes were distributed according to a normal (i.e., Gaussian) distribution centered around a mean of three, whereas the outgoing nodes displayed a very heavy long tail toward higher values, as well as a large number of null values (i.e., leaves, that is, genes only receiving and not sending any signals within the GRN). The long tail was mostly due to 4 genes (CHGA, CHGB, STMN2, and HMGB2). The intensity of all outgoing interactions for these four genes is shown in Figure 4c. HMGB2 displayed a specific pattern with quasi-exclusively inhibitory interactions, whereas the other three genes displayed a relatively similar number of positive and negative interactions. It should be noted that all but HMGB2 are part of the defining gene signatures (see Materials and Methods).

One key feature of our inference process is the possibility of decomposing the overall GRN into dynamical subparts, where each edge appears at the time point transition for which it was detected with the strongest intensity by the inference algorithm ( Figure 5; ¹⁸ see Materials and Methods for a formal definition).

The dynamics start with a small number of genes being turned on by the stimulus ( Figure 5a). In our previous studies, the stimulus was explicitly identified as an exogenous medium change/addition. ^{18, 35} Here, it should be understood as a complex influence exerted by the cell environment, whether hormonal or through cell-cell interactions, which induces SCP cell differentiation. It is needed to “set the GRN in motion” and push the cells out of their quasi-steady state.

There was a much larger number of interactions important for the next time point ( Figure 5b), including a very large set of genes being repressed by HMGB2. All known SCP-specific genes were repressed at that point, which was expected. A CHGA-CHGB self-reinforcing loop activates genes from the chromaffin lineage, while a few sympathoblast-specific genes (STMN2, HAND2) are activated.

At the next time point, an RTN1-PENK toggle switch was apparently supported by positive loops on both sides ( Figure 5c). We can see a lot of genes being repressed by CHGA and CHGB genes. This includes the repression of MPZ-and PLP1 SCP-specific genes by CGHB, a chromaffin-specific gene.

The two latest time points showed a very small number of interactions ( Figure 5d and 5e).

It is clear that the two time points where there is a large number of interactions correspond to the times at which one observes a steep variation in the stemness score ( Figure 2). Altogether, this time decomposition offers a clear illustration of a signal propagating through the GRN in the form of waves of gene activation and repression. ³⁵

Finally, we extracted from the full GRN a subset of the genes connected by specific dynamical motifs, that is, either cross-positive self-reinforcing loops or mutually repressive toggle switches ( Figure 6). Among these, only one direct toggle switch can be observed, the one linking PENK to RTN1 genes. All other motifs were positive. One can see a group of genes (COL5A2, PTPRZ1, EDRNB, POSTN), including one SCP-specific gene (POSTN) which reinforce together and are repressed concomitantly by chromaffin specific genes (CHGA and CHGB) and a sympathoblast-specific gene (STMN2).

One of the main advantages of our approach is that it produces executable GRN models, ³⁶ the time-dependent evolution of which can be simulated, and the resulting simulated dataset can be compared with the experimentally observed dataset.

There are many ways in which the results of a simulation can be assessed. Using a UMAP representation allows for quick evaluation of the quality of the fit. This was first used to calibrate the process times t1–t5 in real time ( Figure 7).

We first tested an evenly spaced time frame, where each time point was set 4 h apart ( Figure 7b). All the cells were instantly projected from their t1 position to the final position of the process. Therefore, we attempted different adjustments and obtained the best visual fit using a time sequence extending up to 80 h ( Figure 7c), which was used for further simulations.

We then assessed the extent to which the simulation captured the time-dependent evolution of a few genes whose expression was characteristic of the three cell types we were trying to reproduce ( Figure 8).

Overall, the gene expression pattern, when observed through its UMAP projection, was quite well reproduced by our GRN simulation. Most notably, chromaffin/sympathoblast bifurcation was quite well captured by our model.

Another way to assess the quality of the fit is to investigate the distance between the distribution of gene expression values of the experimentally observed mRNA distribution and the simulated distribution. As shown in Figure 9, there was a close proximity between the two distributions, as assessed by a low Kantorovich distance for the four genes displayed. Altogether, although some refinement could be performed (see discussion), we considered the fit quality of the experimental dataset to be quite satisfactory.

We finally analyzed the run-to-run variability of the model, given the intrinsically stochastic behavior of our gene expression model (see “GRN simulation” upper). We therefore performed 10 runs of the GRN without specifying the seed, and analyzed the result by taking advantage of the clear separation of the three cell types in the 2D UMAP space ( Figure 1c). This allowed us to count the number of cells belonging to each cell type following the simulation ( Figure 10). We observed a very reproducible behavior of our model with only marginal differences between the number of cells for each cell type between different simulations.