Modeling gene expression cascades during cell state transitions

Summary During cellular processes such as differentiation or response to external stimuli, cells exhibit dynamic changes in their gene expression profiles. Single-cell RNA sequencing (scRNA-seq) can be used to investigate these dynamic changes. To this end, cells are typically ordered along a pseudotemporal trajectory which recapitulates the progression of cells as they transition from one cell state to another. We infer transcriptional dynamics by modeling the gene expression profiles in pseudotemporally ordered cells using a Bayesian inference approach. This enables ordering genes along transcriptional cascades, estimating differences in the timing of gene expression dynamics, and deducing regulatory gene interactions. Here, we apply this approach to scRNA-seq datasets derived from mouse embryonic forebrain and pancreas samples. This analysis demonstrates the utility of the method to derive the ordering of gene dynamics and regulatory relationships critical for proper cellular differentiation and maturation across a variety of developmental contexts.


Highlights
Fitting pseudotimeordered expression profiles to interpretable functional forms

Derivation of transcriptional cascades to define a pseudotime trajectory
Inference of directionality of regulatory interactions

INTRODUCTION
Changes in gene expression underlie the intrinsic molecular processes governing differentiation, enabling cells to change their morphology and function.These changes can occur in part due to extrinsic cues from signaling molecules 1 or temperature and oxygen levels in the organism's environment, 2,3 as well as intrinsic mechanisms such as the asymmetric distribution of cellular components during cell division. 4hese processes result in modifying the expression levels of genes that are critical for cell fate specification, most importantly transcription factors, which can initiate or block the expression of downstream target genes, including other transcription factors.The sequential activation and repression of transcription factors and their target genes can give rise to a cascade of gene expression, whereby an initiating event can regulate a hierarchy of downstream genes essential for the cell to acquire subsequent cell states.For example, the Pax6 / Eomes / Tbr1 transcription factor cascade directs the progression of radial glia to intermediate progenitor to postmitotic projection neuron in the developing cortex, 5,6 and the transcription factor cascade initiated by Neurog3 controls the differentiation of endocrine progenitor cells to mature pancreatic cells. 7,8It is therefore critical to accurately deduce gene expression cascades in order to determine which genes are responsible for specific cell fate changes during differentiation and maturation.Single-cell RNA sequencing (scRNA-seq) enables sampling the gene expression profile of thousands of cells in an individual sample.However, it is necessary to destroy the cell in order to measure its transcriptome, thereby making it impossible to observe how the cell and its gene expression profile would have altered in the future.Nonetheless, it is possible to order cells along a trajectory which accurately recapitulates the progression of cells as they transition from one cell state to another.This ordering of cells along a trajectory is known as pseudotime, which is essentially a mapping of single-cell transcriptomes to a developmental timeline.][11][12][13][14] Based on the ordering of cells along a pseudotemporal trajectory, it is possible to measure the dynamics of gene expression as cells undergo cell state transitions.Current algorithms typically model gene expression dynamics along pseudotemporal trajectories by fitting their expression profiles using generalized linear models, 12,15,16 with the ultimate goal of determining if gene expression significantly varies as a function of pseudotime.Other methods attempt to deduce pseudotime-dependent gene interactions by calculating a similarity measure between the expression levels of the ''present'' of one gene, and the ''past'' of another gene using correlation 17 or mutual information. 18However, these methods do not calculate an explicit ordering of expression dynamics along a pseudotime trajectory, and require user-defined cutoffs for determining meaningful interactions.
Here, we present a method to better understand the cascade of gene expression dynamics underlying cell state transitions.We are interested in answering questions such as, if two genes are up-regulated during a cell state transition, is one gene up-regulated before the other, or are they up-regulated simultaneously?Furthermore, is it possible to estimate a certainty in the timing of their expression dynamics?In this paper, we address these questions by explicitly modeling gene expression over a pseudotime trajectory using a set of functions that reflect biological state switches, and that model the dynamic behaviors of gene expression within cells as they differentiate.We formulate the problem using a Bayesian inference framework and use an ensemble sampler Monte Carlo Markov chain (MCMC) approach 19 to sample from the posterior distributions over the parameter spaces of the various functions, and determine which model best fits the data.This provides an explicit ordering of genes along a pseudotemporal trajectory based on inflection point estimates, enabling the description of expression dynamics in terms of transcriptional cascades, estimating differences in switch times of gene expression, and annotation of potentially causal gene interactions in gene regulatory networks.
We will introduce our modeling framework in general terms in the first section of the results.A more detailed description is provided in the STAR Methods section.We then apply our method in multiple developmental settings, in which we dissect the transcription factor cascades underlying cortical neurogenesis and pancreatic beta cell development across multiple scRNA-seq datasets.We also show how our method can be used to infer potential upstream regulators of a given gene of interest.Finally, we utilize our method to deduce the gene expression cascade of the Notch signaling pathway in the developing cortex in order to highlight the applicability of our method to gene sets beyond transcription factors.These examples demonstrate the ability of our method to accurately model the dynamics of gene expression during cell state transitions, and highlight the biological insights our method enables.

Modeling gene expression dynamics along pseudotime trajectories
The goal of the method presented here is to decide if a state switch (up-to down-regulation or down-to up-regulation) occurs along a pseudotemporal trajectory, and at what pseudotime these switches occur, in order to determine the timing and ordering of activation and repression during cell state transitions.In order to do this, we first define a set of functions which can model a wide variety of expression dynamics, and for which state changes are well defined and interpretable, namely at the inflection points of each function.The functions are then fit to the normalized expression levels for each gene across cells ordered by their relative pseudotempoal ordering.The functions used for fitting are defined as follows, (Equation 1) Here, f unif is a uniform function with b > 0, which models the absence of dynamics in gene expression along a pseudotime trajectory.f gauss is a Gaussian function with parameter constraints a > 0; b > 0, s > 0; and 1 % t 0 % N, with N = number of cells in the pseudotime trajectory.f sig is a sigmoidal function with parameter constraints L > 0; b > 0; and 1 % t 0 % N. Finally, f dsig is a double sigmoidal function with the formulation described in the study by Baione et al. 20 and parameter constraints b min > 0; b mid > 0; b max > 0; k 1 > 0; k 2 > 0; and 1 % t 1 < t 2 % N. The motivation for using these functions is based on observations from biological scenarios during development. 21For instance, during differentiation, genes can display a shift from one steady state to another, which can be modeled using a sigmoidal function.They can also exhibit impulse patterns of up-regulation followed by a return to basal levels, which can be modeled using a Gaussian function.Finally, double sigmoidal functions can model impulse patterns with asymmetric increase and decrease rates and different initial and terminal basal levels, as well as stepwise up and stepwise down expression patterns (Figure S1).We formulate the problem of fitting gene expression profiles in cells ordered along a pseudotime trajectory as a Bayesian inference problem, and estimate parameters for each function using an ensemble sampler MCMC approach 19 (see STAR Methods).Based on the best-fitting function to the gene expression profiles, genes are ordered according to the relative occurrence of inflection point estimates to provide temporal estimates of gene expression cascades, and regulatory interactions between genes are deduced, enabling a detailed characterization of the molecular processes underlying cellular transitions.

Transcriptional cascades during cortical neuron differentiation
We first applied our method to differentiating forebrain dorsal neural stem cells during mouse development at embryonic stage e13.5.The input to the method consists of a set of cells ordered by pseudotime, t = 1; .;N,and the expression levels (counts) of genes within those cells.Cells from the Atlas of the Developing Mouse Brain 22 were initially subset to non-dividing forebrain dorsal cells consisting of neural stem cells, intermediate progenitors (IPs), and neurons at embryonic stage e13.5.A pseudotime ordering was estimated using diffusion pseudotime 9 (Figure S2).All dividing cells were excluded for the pseudotime estimation due to their expression of a transcriptional program that is independent of the underlying cell type, potentially confounding pseudotime estimates.
In differentiating cells along the mouse e13.5 forebrain dorsal neural stem cell (NSC) / IP / neuron trajectory, 60 out of 510 (11.8%) transcription factors (derived from the study by Lambert et al. 23 ) that were expressed in at least 1% of cells had a non-uniform fit (Figure 1; Table S1).Initially, Gli3, a gene that is required for maintaining cortical progenitors in active cell cycle, 24 was down-regulated in a state-switch manner with a sigmoidal fit, along with Sox9 and Hes1, which are both required for neural stem cell maintenance. 25,26Subsequently, other genes important for neural stem cell maintenance including Sox1, Sox2, Hes5, and Pax6 were down-regulated.Genes exhibiting a state-switch or stepwise up-regulation included Neurod2, Sox11, and Neurod6, which play a critical role in inducing cell-cycle arrest and neurogenic differentiation in the developing cortex, [27][28][29] followed by Tbr1 and Bcl11b, markers of deep-layer cortical neurons generated during early cortical neurogenesis.Subsequently, Satb2 and Bhlhe22, markers of upper-layer cortical neurons generated during later stages of neurogenesis, 30 were up-regulated.Interestingly, four transcription factors were found to be transiently down-regulated using a double sigmoidal fit, including Mycn, Jun, Ybx1, and Jund.Genes exhibiting a transient up-regulation (Gaussian or double sigmoidal fit) included Hes6 and Eomes, markers of cortical IPs, 31 as well as Neurog2 and Sox4, which are required for IP cell specification and maintenance via activation of Eomes. 32hese results demonstrate that the functions which best fit the expression profiles of dynamically expressed genes (genes exhibiting a nonuniform fit) largely reflect the known biological role these genes play during differentiation.Furthermore, the relative ordering of inflection point estimates for dynamically expressed transcription factors along the mouse e13.5 forebrain dorsal NSC / IP / neuron trajectory accurately recapitulates known temporal orderings that are essential for the differentiation of cortical neurons.Finally, in order to justify the functional forms we used, we performed a PCA of the gene expression profiles.Genes with a non-uniform fit fill the extremes of the principal component space (Figure S3), indicating that the functional forms we used to model the pseudotime-ordered gene expression profiles are able to capture most of the variability in the data.

Constructing regulatory interactions during cortical neurogenesis
We then compared a set of transcription factors forming an essential regulatory network underlying cortical neuronal differentiation including Pax6, Neurog2, Eomes, and Tbr1, 33 as well as the neural lineage bHLH factor, Neurod4 (Figure 2A).Neurog2 and Eomes exhibited a transient up-regulation, with both genes having a double sigmoidal fit.Pax6 and Tbr1 were fit using a sigmoidal function, with Pax6 exhibiting a stateswitch from high to low expression, and Tbr1 from low to high expression.Neurod4 was fit using a Gaussian function, and was specifically expressed transiently in mid-stage Eomes + cells.These genes were then ordered according to the pseudotemporal occurrence of inflection point estimates (Figure 2B), whereby Neurog2 was found to be up-regulated before Eomes, followed by the up-regulation of Neurod4 and down-regulation of Pax6.Subsequently, Tbr1 was up-regulated, followed by down-regulation of Neurod4, Neurog2, and finally Eomes.Neurod4 exhibited a brief, transient impulse expression pattern within mid-stage Eomes + cells, reflecting previously studied expression patterns of Neurod4, which is only expressed in a subset of Eomes + cells in the mouse e14.5 cortex. 34y comparing inflection point estimates of these genes (see STAR Methods), we were able to reconstruct previously validated regulatory interactions (Figure 2C).The initial up-regulation of Neurog2 just before Eomes up-regulation suggests that Neurog2 initiates expression of Eomes in intermediate progenitors.This relationship has been shown in mouse e13 embryos via electroporation of Neurog2 cDNA into the ganglionic eminence, where both Neurog2 and Eomes are not expressed, resulting in ectopic expression of Eomes. 35Neurog2 has also been shown to directly activate Neurod4 in cortical IP cells using a luciferase reporter assay, 36 which we also recapitulate based on the sequential up-regulation of Neurog2 and Neurod4.Furthemore, it has been shown that both Neurog2 and Eomes induce Tbr1 expression, 36 which we also infer based on the up-regulation of Tbr1 following both Neurog2 and Eomes.Interestingly, directly after Eomes and Neurog2 were upregulated, Pax6 was down-regulated, suggesting a negative feedback loop, whereby Pax6 activates both Eomes and Neurog2, which then both in turn repress Pax6, a relationship which has been previously described in the developing mouse cortex. 37

Inferring shared upstream regulators of Eomes
We next explored potential upstream regulators of Eomes in mouse e13.5 forebrain dorsal cells across two samples in order to deduce high confidence regulators of Eomes and determine how robust our method is across biological replicates.We applied our method to forebrain dorsal cells in a mouse e13.5 biological replicate (Figure S4; Table S2).Transcription factors with a positive inflection point occurring simultaneously with or before the first inflection point of Eomes, as well as those with a negative inflection point occurring after the first inflection point of Eomes, were labeled as positive upstream regulators.We furthermore included all co-activators and co-repressors (derived from the study by Siddappa et al. 38 ) that exhibited a transient up-regulation, with the first inflection point occurring simultaneously with or before the first inflection point of Eomes.In total, 25 positive upstream regulators were found in the first sample, and 27 were found in the second sample, with an overlap of 21 genes across the two (Figure 3A; Figure S5).Furthermore, the relative ordering of inflection points of these genes along the cortical differentiation trajectory strongly agrees across both datasets, with one exception being Tfap2c, which was fit to a sigmoidal function in the first sample, and Gaussian function in the second sample.
Within the set of inferred transcription factors regulating Eomes expression were Neurog2 and Pax6, which are known to directly activate Eomes in the developing mouse neocortex, as described in the previous section.The co-regulators Dll1, a key ligand for activating Notch signaling, and Chd7, a chromatin remodeler, have also been implicated in the formation of IP cells, 39,40 although their role as a co-activator of Eomes has not been established to our knowledge.These results validate the utility of our method in discovering upstream regulators of a given gene of interest.The remaining potential activators of Eomes warrant further experimental validation.
Furthermore, the genes that repress Eomes in maturing IP cells, thereby enabling the differentiation of these cell types into neurons, are largely unknown. 33The transcription factor Mycn, a gene critical for normal brain development, 41 has been shown to down-regulate Eomes in neuroblastoma cell lines 42 ; however, its role in regulating Eomes expression in maturing IP cells is not well understood.In differentiating cells along the forebrain dorsal NSC / IP / neuron trajectory in both mouse e13.5 samples, Mycn was expressed in a transient down-regulation pattern and best fit using a double sigmoidal function (Figures 3B and 3C).In both samples, Mycn up-regulation occurred simultaneously with Eomes down-regulation, signifying that Mycn may play a role in the differentiation of cortical neurons by down-regulating Eomes in maturing IPs.

Dissecting Notch signaling during cortical neurogenesis
To demonstrate the applicability of our method to genes beyond transcription factors, we investigated the dynamics of Notch signaling along the forebrain dorsal NSC / IP / neuron trajectory in e13.5 mouse embryos.Shared dynamically expressed genes involving ligand-receptor pairs of Notch receptors from the study by Shao et al. 43 in both embryonic samples were estimated (Figure 4).In both samples, Mfap2, which can interact with the extracellular domain of Notch1, 44 however whose role is poorly understood in the regulation and differentiation of cortical NSCs, was up-regulated within forebrain dorsal NSCs, and down-regulated in neuronal cells.This indicates that Mfap2 may play a general role in Notch signaling within differentiating cortical NSCs, whose actions are not specific to a given cell type.Dll1 was up-regulated in early IPs, followed by the up-regulation of Dll3 in later stage IPs, confirming the selective basal expression of Dll3 from in vivo studies. 33Furthermore, Mfng, a glycosyltransferase which increases the ability of Notch1 to bind to Dll1, 45 was up-regulated shortly after Dll1 up-regulation in both samples within IPs, indicating that this gene becomes activated sequentially after the activation of Dll1.Dll1 was then down-regulated within IPs, suggesting that this gene is not essential for further IP differentiation into neurons.Finally, Notch1 was down-regulated in maturing IPs, followed by down-regulation of Mfap2, Mfng, and Dll3 in neurons.These results highlight the ability of our method to dissect the complex dynamics of signaling pathways within differentiating cell types.

Transcriptional cascades in mouse pancreatic beta cell development
To demonstrate the utility of our method in other developmental contexts, we applied our method to a scRNA-seq dataset of pancreatic cells derived from mouse e14.5 embryos, 46 subsetting to cells belonging to the beta cell lineage.When measuring the expression dynamics of a set of genes known to play an essential role in the specification and maturation of pancreatic beta cells, 8 we find a well-defined transcriptional cascade which largely agrees with previously characterized gene expression cascades (Figure 5A).Interestingly, we find one exception to this cascade, Neurod1, which is up-regulated at a later stage of beta cell maturation than previously reported (Figures 5B and 5C).We are also able to measure the sequential up-regulation of Pax6 and Pdx1, followed by Mnx1, and ending with the insulin gene expression regulator Isl1, thereby providing a more explicit ordering of the expression cascade in maturing beta cells than previously established.Furthermore, with this approach, we can model the expression dynamics of all transcription factors (Figure S6; Table S3), enabling a detailed overview of the full gene expression cascade underlying pancreatic beta cell differentiation.

DISCUSSION
In this paper, we explored an approach to model the gene expression dynamics in cells ordered by a pseudotime trajectory using a fully Bayesian framework.This framework enabled us to fit the gene expression profiles of cells undergoing cell state transitions to a set of functions that are able to model complex transcriptional dynamics.From these fits, we were able to order genes along a gene expression cascade which describes the molecular dynamics underlying cell state transitions, and deduce regulatory interactions.
We first applied the method to differentiating forebrain dorsal neural stem cells into neurons in mouse e13.5 embryos.By ordering transcription factors by the relative occurrence of inflection point estimates, we were able to reconstruct the transcriptional cascades underlying neuronal differentiation within the developing cortex, and model the dynamics of gene expression for all genes along the trajectory.However, genes can undergo further dynamic changes including post-transcriptional and post-translational modifications, and localization changes within the cell, all of which can have a large impact on function and regulation.While transcriptomics data are unable to identify these changes, the dynamics we uncover from gene expression data can still shed light on their regulatory roles.
By comparing the relative timing of expression dynamics of the transcription factors Pax6, Neurog2, Eomes, Neurod4, and Tbr1, which form a regulatory network underlying cortical neuron differentiation, we were able to infer known causal interactions.However, reconstructing a gene regulatory network using all genes with a non-uniform fit would lead to many false positives, in part due to the simultaneous activation of multiple pathways involving different genes.Thus, we believe one of the main utilities of our approach is to infer the directionality of regulatory interactions, especially in cases where an interaction has been measured but the directionality is unknown.
We then identified potential upstream positive regulators of Eomes, an essential gene for the formation of IPs.Subsetting to genes which have similar dynamics across biological replicates revealed a set of high-confidence potential upstream regulators.Not only did we recover validated activators of Eomes, such as Pax6 and Neurog2, but we also detected a number of other transcription factors whose roles in Eomes activation have not been fully characterized.The enrichment of known DNA-binding motifs of these transcription factors in the promoter and enhancer regions of Eomes may provide further evidence for the regulatory role of these genes in Eomes expression.We also identified a potential negative regulator of Eomes, the transcription factor Mycn, whose role in cortical IP maturation has not been fully explored.Wet lab experiments, such as knockin or knockout experiments, or chromatin immunoprecipitation sequencing experiments, would need to be performed in order to validate the roles of these transcription factors in the regulation of Eomes expression.
We further demonstrated the applicability of our method to genes beyond transcription factors by comparing the expression dynamics of genes involved in the Notch signaling pathway.This analysis revealed a sequential up-regulation of the Notch receptor ligand Dll1 in early IPs, followed by Mfng, and finally Dll3 in maturing IPs.This activation cascade supported the selective expression of Dll1 and Dll3 in apical and basal IPs, respectively, further demonstrating the utility of comparing genes according to inflection point estimates to dissect signaling pathways.
We also applied our method to differentiating pancreatic beta cells in mouse e14.5 embryos.Based on this analysis, we were able to reconstruct a gene expression cascade that defines beta cell maturation.In this analysis, we highlighted a gene that deviated from the established literature, Neurod1, whose up-regulation along the cascade occurred later during beta cell development than previously established.Followup experiments are needed to validate these findings.
In order to place our method in a broader context, we compared our results with Monocle 3 12 and tradeSeq, 16 which perform statistical tests to determine if a gene is differentially expressed along a pseudotime trajectory, in cells from the e13.5 forebrain dorsal NSC / IP / neuron trajectory.While the overwhelming majority of genes with a non-uniform fit from our method were also found to be significantly differentially expressed by these two methods, both methods detected at least six times more genes to be significant compared to our method (Figure S7).Thus, we conclude that our method is more stringent in detecting genes exhibiting dynamic changes along a trajectory.Furthermore, while the relative ordering of gene expression dynamics along a trajectory is not readily available using these two methods, we are able to explicitly infer this using our method based on inflection point estimates.Similar to our method, the authors of the original diffusion pseudotime publication used derivative estimates of smoothed gene expression profiles to order gene dynamics along a pseudotime trajectory. 9owever, the authors only used derivative estimates to measure switch-like transitions and not transient up or down transitions, and only provide point estimates of these transitions.We are able to model a higher variety of transitions, and based on the MCMC samplings, quantify the uncertainty in the timing of these transitions using the posterior distribution of the parameter fits.
To measure the dependence of our method on the pseudotime method used to order cells, we ran our method on the pseudotime-ordered cells from the e13.5 forebrain dorsal NSC / IP / neuron trajectory using both Slingshot 10 and Monocle 3, 12 and compared them with the diffusion pseudotime estimates (Figure S8).Overall, the fits were largely consistent independent of the pseudotime method used to order the cells, indicating that our method is robust to fluctuations in pseudotime estimates and underlying pseudotime method.
While we focused specifically on cells along the forebrain dorsal NSC / IP / neuron trajectory, and pancreatic beta cell development, the method presented in this paper can be applied to any scRNA-seq dataset where cells can be ordered along a pseudotime trajectory.Our method is able to reconstruct transcriptional cascades in order to deduce critical genes for cell state transitions.It is also able to predict regulatory interactions, as well as gene interactions involved in different signaling pathways.Therefore, we believe this approach can provide useful insights into the molecular underpinnings involved in a variety of developmental biology contexts.exhibiting high expression levels of G2M cell cycle genes were subsequently filtered, as well as clusters with a subpallium (ventral cortical) identity, hippocampal identity, and Cajal-Retzius neurons.The above procedure was re-run until the only subsequent populations in the sample consisted of forebrain dorsal NSCs, IP cells, or neurons based on the expression of known marker genes for the respective populations.Diffusion pseudotime estimates 9 for each cell were then estimated was after running a diffusion map embedding and assigning a starting cell.The raw count data across all cells ordered by diffusion pseudotime were then stored and the MCMC procedure was run on the resulting count matrix.

pancreas development samples
The raw count data for the pancreas endocrinogenesis dataset 46 was downloaded from http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132188.The raw count data was loaded into scanpy 47 for downstream analyses.Cells were then initially subset to samples corresponding to e14.5 embryos.All cells with a positive G2M score in the metadata were initially filtered.Following this, the count data was processed in a similar fashion to the Atlas of the Developing Mouse Brain dataset using scanpy.Diffusion pseudotime estimates were calculated and the raw count data across all cells ordered by diffusion pseudotime were then stored and the MCMC procedure was run on the resulting count matrix.

Establishing a likelihood model
The negative binomial distribution has been shown to accurately describe the count data generated in scRNA-Seq experiments without the need to account for zero-inflation resulting from ''dropout'' events. 51The probability mass function for the negative binomial distribution can be parameterized using the mean, m ˛R+ , and dispersion parameter, 4 ˛R+ , with y ˛N, as follows, (Equation 2) The mean and variance of the random variable Y $ NBðm; 4Þ which follows a negative binomial distribution is then E½Y = m and Var½Y = m + m 2 4 .For a gene g with measured counts of Y !g = fy gt g t = 1;::;N along a pseudotime trajectory with fixed pseudotime-step interval, m !g = fm gt g t = 1;.;Nand 4 !g = f4 gt g t = 1;.;N the mean and dispersion at corresponding pseudotimes, the full likelihood of observing Y where pðy gt m gt ; 4 gt Þ is the negative binomial probability mass function.The full log-likelihood is then: It was shown that when fitting scRNA-Seq UMI count data to a negative binomial model, data are consistent with a global dispersion parameter independent of the expression level of a given gene, and that fitting a dispersion parameter to each gene individually leads to overfitting. 52Therefore, a global estimate of 4 can be used for every gene independent of pseudotime, and 4 g != f4 gt g t = 1;.;N is replaced with a constant 4 in Equation 4. A dataset specific 4 using genes which exhibit lower levels of overdispersion is estimated, since the expression levels in these genes reflect the technical rather than the biological variability.To do this, the log10 mean counts for each gene are binned into five equally spaced bins, and a linear fit between log10 mean and log10 variance of counts in each bin is estimated.Genes within the top 20th percentile of the difference between the estimated variance and the expected variance using the linear fit in each bin are then filtered.The remaining genes are used to fit the non-linear relationship between the mean (m) and variance (s 2 = m + m 2 4 ) using unconstrained non-linear least squares (Figure S9).
Here, 4 estimates the dispersion based on genes which do not exhibit high variability in the dataset, and therefore captures the technical variability in the dataset.This technical variability is in large part driven by the varying number of UMI counts captured in each cell, as well as other factors including library quality and amplification bias.Thus, the full log-likelihood of observing counts Y !g = fy gt g t = 1;::;N for gene g along a pseudotime trajectory given the mean at corresponding pseudotime points m !g = fm gt g t = 1;.;N, becomes, 5) where 4 is a global parameter estimated using the procedure described above.For scRNA-Seq methods which sequence only from one end of the transcript and not full-length protocols, normalization does not need to account for the total transcript length.In this case, for a given cell i, let M i be the number of UMIs in cell i, and y gi be the number of UMIs for gene g in cell i.In this paper, we use the median number of UMIs across all cells in the dataset as a size factor M, that is, M = medfM i g i = 1;.;N .
Then, the log-normalized expression levels for gene g in cell i is defined by the following mapping, h y gi = ỹgi = ln y gi M i M + 1 : (Equation 6) The functions ðf unif ; f gauss ; f sig ; f dsig Þ described in Equation 1 are then fit to the pseudotemporally ordered expression profile for gene g, fỹ gt g t = 1;.;N, in the log-normalized expression space with the objective function to maximize defined by the likelihood in Equation 5.The means m !g = fm gt g t = 1;.;Nare then calculated by mapping the function values evaluated at t = 1; .;N back to count space using the inverse of Equation 6.The full log-likelihood estimate is then evaluated by plugging in the m !g values and global estimate for 4 into Equation 5.This procedure can be summarized as follows.We want to solve for f a ðt; qÞ, which maximizes the following likelihood,  7) where f a ˛ðf unif ; f gauss ; f sig ; f dsig Þ.

Model inference using MCMC
Under the framework presented above, solving for f a ðt; qÞÞ can be formulated as a Bayesian inference problem, which we solve using an ensemble sampler MCMC approach. 19This provides an estimate of the posterior distribution over the parameter space for each of the parameters in the different functions ðf unif ; f gauss ; f sig ; f dsig Þ described in Equation 1.For each of the models, the priors used for the different parameters are summarized in Table S4.Note, in Table S4 (Equation 8) The uniform priors in Table S4 are uninformative, however, they provide bounds on the parameters to keep them in interpretable and meaningful ranges.The slope parameters k in the sigmoidal function, and k 1 and k 2 in the double sigmoidal function, have a folded normal prior with 0-mean and 0.1 variance, which is used to ensure that the slope has a low magnitude.This prior is used because differences in the function once the slope becomes relatively large are minimal.Finally, the folded normal prior on s in the Gaussian with 0-mean and N= 10 variance is used to ensure that the curve does not become very flat.
In this paper, we use the ensemble sampler MCMC proposed by Goodman & Weare in 2010 19 with implementation by Foreman-Mackey et al. 53 An initial guess is needed as a starting point from which a walker begins in the ensemble sampler.For the Gaussian and sigmoidal functions, initial guesses are derived from a non-linear least squares fit for each function on the log-normalized pseudotime expression levels using scipy's 'curve_fit' function, with added Gaussian noise.For the double sigmoidal function, initial guesses are randomly chosen to cover the varieties of different forms the functions can have.For the uniform function, initial guesses are randomly chosen from a uniform distribution over the interval 0.01 and maximum expression level for the gene of interest.The number of walkers used is four times the number of parameters for each function -28 for the double sigmoidal fit, 16 for the Gaussian fit, 16 for the sigmoidal fit, and 4 for the uniform fit.This enables a wide sampling across the search space of parameters.
The MCMC is then run for a total of 10,000 iterations.There is generally no consensus on how many iterations to run an MCMC algorithm. 53housands of iterations are typically desirable to allow the process to reach a steady-state.After reaching the steady-state, the MCMC will sample from the posterior distribution over the parameter space, enabling an estimate of the posterior distribution for each parameter.Iterations before reaching the steady-state are discarded, as these are not sampled from the target distribution.This is called the ''burn-in'' phase.For this implementation, a burn-in of 5; 000 iterations was used (Figure S10).Some MCMC walkers can get stuck near a local maximum.These walkers typically have a low acceptance rate, that is the proportion of moves for which the MCMC sampler generated parameter values that differed from the previous sample.One common practice is to prune these walkers from the final MCMC output.For example, walkers can be pruned which get stuck in irrelevant local optima by clustering the likelihood of the walkers and removing the clusters with lower likelihoods. 54For this implementation, half of the MCMC walkers are pruned with the lowest acceptance rate in order to remove potentially stuck walkers (Figure S11).

Model selection
We use a probabilistic model selection technique, the bayesian information criterion (BIC) 55 to score the different models, and select the model with the best score.The BIC is defined as follows, BIC = k lnðnÞ À 2 lnð b LÞ; (Equation 9) where n = number of data points, k = number of parameters in the model, and b L = maximized value of the likelihood function.In the original formulation of the BIC, the value b L was derived from maximum likelihood estimation.When using an MCMC for model inference, the output consists of a sampling or distribution over the parameter space.It is advantageous to use a likelihood estimate which more closely reflects the optimal parameter regime estimated from the MCMC instead of the parameter regime which maximizes the likelihood.To this end, b L in the Here, T ˛½0; 1000 enables an accurate estimate of b t f under the assumption that r f ðtÞ approaches 0 by t = T for each parameter.The autocorrelation function (Figure S13) and autocorrelation time (Figure S14) is estimated for each parameter separately.
For a general comparison, the autocorrelation times were estimated for all genes using the model with the best fit in the mouse e13.5 forebrain sample (Figure S15).The autocorrelation times increase with the complexity of the model (i.e.number of parameters specified in each model).This is in part expected, since a model with more parameters will generally have a lower acceptance rate due to the higher number of dimensions in which the MCMC has to make proposal moves, leading to higher autocorrelations for each parameter.Nonetheless, the autocorrelation times are fairly robust for each model.
Thinning is an approach to use every k-th iteration of the MCMC walkers, where k = t f would represent an i.i.d.sampling of the posterior distribution.However, various publications indicate that thinning is often unnecessary and results in reduced precision. 57,58Therefore, no thinning of the MCMC walkers was used in this analysis.
Another way to visualize the posterior distribution over the parameter space derived from an MCMC is a corner plot (Figure S16).The corner plot highlights both the two dimensional projections over the parameter space across iterations of the MCMC, as well as the marginal posterior distribution for each individual parameter (highlighted in the upper plots).Some parameters are more correlated with each other than others, indicating underlying covariates within the model parameters.However, the marginal posterior distributions do not appear to be multimodal.
These heuristics provide some insight into the ability of the ensemble MCMC sampler to provide an accurate sampling of the posterior distribution over the parameter space.

Estimating inflection points
Inflection points occur where the curvature of a function changes sign.At inflection points, the first-order derivative, or rate of change, of a function reaches a local maximum or local minimum.At an inflection point, the second-derivative of a function passes through 0 with the second derivative changing sign from positive (concave upward) to negative (concave downward) or vice versa.The inflection points of the Gaussian, sigmoidal and double sigmoidal fits can be used to compare the relative timing of when genes exhibit a state transition along a pseudotime trajectory.To estimate the inflection points of the different functions, first solve for x at which the second-derivative of the function is zero.For the Gaussian function, f gauss ðtÞ, sigmoidal function f sig ðtÞ, and double simgoidal function f dsig ðtÞ defined in Equation 1, the second derivatives are f 00 gauss ðtÞ = a s 4 e À ðt À t 0 Þ 2 2s 2 ðt À ðt 0 À sÞÞðt À ðt 0 + sÞÞ; For the Gaussian function, f gauss ðtÞ, two inflection points occur at t ˛ðt 0 À s;t 0 + sÞ.For the sigmoidal function, f sig ðtÞ, one inflection point occurs at t = t 0 .The estimates for the inflection points are then measured from the parameters ðt 0 À s; t 0 + sÞ for the case of the Gaussian and t 0 for the case of sigmoidal function at each MCMC iteration.Finally, for the double sigmoidal function, f dsig ðtÞ, the number of inflection points can vary.However, if all parameters are fixed besides k 1 , then, f 00 dsig ðtÞ/0 as k 1 increases.Similarly, if all parameters are fixed besides t 1 , then f 00 dsig ðtÞ/0 as t 1 decreases.That is, for k 1 [ 0, i.e. the transition from b min to b mid occurs rapidly, then an inflection point will occur very close to t 1 .Similarly, for k 2 [ 0, i.e. the transition from b mid to b max occurs rapidly, then an inflection point will occur very close to t 2 .Also, the further apart t 1 and t 2 are from each other, the closer the inflection points are to t 1 and t 2 .To ensure the inflection points occur very close to t 1 and t 2 , at each iteration of the MCMC, a move is only accepted in cases where signðf 00 dsig ðt 1 À dtÞÞ Ã signðf 00 dsig ðt 1 + dtÞÞ < 0 and signðf 00 dsig ðt 2 À dtÞÞ Ã signðf 00 dsig ðt 2 + dtÞÞ < 0 for dt = 1.The estimates for the inflection points are then calculated from the parameters t 1 and t 2 at each MCMC iteration.

Comparing inflection points
Regulatory interactions were inferred based on the relative timing of inflection point estimates (Figure 2).If there was an overlap of at least 1% in the inflection point estimates between two genes across MCMC iterations, then these were assumed to have a simultaneous switch state.A regulatory interaction between the two was mutually positive if the inflection points had the same sign, and mutually negative if the inflection points differed in sign.The overlap between two inflection points is estimated by binning the inflection point estimates across all MCMC iterations to 100 equally spaced bins, starting at the minimum inflection point estimate across both genes and ending at the maximum inflection point estimate across both genes.Let fx i g i ˛½1;100 represent this binning domain.If p A ðx i Þ is the percent of counts in the histogram in bin x i for gene A, and p B ðx i Þ is the percent of counts in the histogram in bin x i for gene B, then the overlap between the two, PðA = BÞ, is

Figure 1 .
Figure 1.Transcriptional cascades in mouse 13.5 forebrain dorsal cells (A) Gene expression profiles of transcription factors with non-uniform fits are displayed as a heatmap.Genes are grouped according to a state-switch from high to low expression (sigmoidal fit) or stepwise down-regulation (double sigmoidal fit), a state-switch from low to high expression (sigmoidal fit) or stepwise upregulation (double sigmoidal fit), a transient up (Gaussian or double sigmoidal fit) expression pattern, and transient down (double sigmoidal fit) expression pattern.(B) The inflection point estimates are shown for the same genes as in (A).Inflection point estimates from double sigmoidal fits are shown in light blue and light red, and those from Gaussian and sigmoidal fits in blue and red.

Figure 2 .
Figure 2. Reconstructing regulatory interactions during mouse e13.5 cortical development (A) Normalized expression levels of essential genes -Pax6, Neurog2, Eomes, and Tbr1 -forming a regulatory network underlying cortical neuron differentiation, as well as the neural lineage bHLH factor, Neurod4, across pseudotime-ordered cells are shown.The curves display a random sampling of the parameters from 100 iterations of the MCMC traces for the best-fitting model for each gene.(B) Inflection point estimates for the genes highlighted in (A).(C) A reconstructed gene regulatory network based on the comparison of inflection points.Positive regulatory interactions which have previously been validated are highlighted as a green solid line, and those which have not been validated as a green dashed line.Similarly, negative regulatory interactions which have previously been validated are highlighted as a red solid line, and those which have not been validated as a red dashed line.

Figure 3 .
Figure 3. Inferring upstream regulators of Eomes across mouse e13.5 embryos (A) The left and right plots show a transcriptional cascade of the shared potential positive regulators of Eomes in forebrain dorsal cells of mouse e13.5 embryos across biological replicates.Transcriptional co-activators and co-repressors (derived from the study by Siddappa et al. 38 ) are shown in orange, and transcription factors (derived from the study by Lambert et al. 23 ) are shown in black.(B) The left panel in the plot displays a random sampling of the parameters from 100 iterations of the MCMC traces for the genes Eomes and Mycn using the double sigmoidal model, the best-fitting model for both genes.The full range of first and second inflection point estimates for both genes is highlighted as a shaded region, with blue indicating a negative inflection point and red a positive inflection point.The middle and right panels highlight the distribution of first and second inflection point estimates across MCMC iterations, respectively.The right panel highlights the distribution of second inflection point estimates across MCMC iterations.p values were estimated as the percentage of overlapping inflection point estimates across both genes after binning the inflection point estimates across all MCMC iterations to 100 equally spaced bins, starting at the minimum inflection point estimate and ending at the maximum inflection point estimate across both genes.(C) The same plot for (B) in cortical cells of the biological replicate.

Figure 4 .
Figure 4. Notch signaling cascade in mouse e13.5 embryos The left and right plots show a transcriptional cascade of the shared ligand-receptor pairs involved in Notch signaling in cells along the forebrain dorsal NSC / IP / neuron trajectories in mouse e13.5 embryos across biological replicates.Annotated cell types are highlighted below.

Figure 5 .
Figure 5. Gene expression cascades in developing mouse e14.5 pancreatic beta cells (A) Schematic diagram of the previously characterized gene expression cascade in developing pancreatic beta cells, based on the study by Wilson et al. 8 (B) The heatmap in the upper panel highlights the expression profiles of transcription factors ordered by the occurrence of their first inflection points.Inflection point estimates are highlighted in the plot below using the same ordering, with double sigmoidal fits shown in light blue and light red, and those from Gaussian and sigmoidal fits in blue and red.The annotated cell type for each cell in the trajectory is highlighted in the middle.(C) Modified gene expression cascade based on inflection point estimates from (B).
, the folded normal distribution is parameterized by m > 0 and s > 0 with probability density function,