Gene Regulatory Network Inference from Single-Cell Data Using Multivariate Information Measures

1Precisely controlled patterns of gene expression are essential for the survival and reproduction of all life forms. Development provides the canonical example, where changes in gene regulation determine the path by which a complex multicellular organism emerges from a single fertilized egg cell. Intricate networks of transcriptional activators and repressors have evolved to regulate the spatial and temporal expression of genes, enabling organisms to adjust transcription levels in response to environmental, developmental, and physiological cues. Elucidating the structure of such gene regulatory networks (GRNs) has been a central goal of much recent systems biology research, and it is now becoming a pivotal stepping stone in dissecting the molecular contributions of complex diseases.

2The structure of GRNs alone does not fully constrain their function, but it serves as an important starting point for further analysis. The simplest mathematical representations of GRNs are static, undirected graphs, where each node represents a gene, and edges depict relationships between transcription factors and their targets. Although by this definition the GRN encapsulates every transcriptional regulatory relationship that could occur within a given organism, this is not a very helpful perspective: clearly, it is not the case that every possible interaction does occur in every cell — for example, the downstream interactions of a transcription factor only occur when it is expressed — hence we can define subsets of the GRN that are active in particular cells and contexts. The structure and dynamics of these active GRN subsets give rise to distinct mRNA expression profiles, and it has been suggested that characteristic expression profiles in different cell types (and under different conditions) result from different stable states of the GRN.

3The introduction of efficient high-throughput expression assays has driven interest in network inference methods that apply statistical approaches to identify likely regulatory relationships between genes based on their expression patterns and potential GRN structures. In addition to correlation-based networks (perhaps the simplest way of identifying putative relationships), Gaussian graphical models, (dynamical) Bayesian networks, regression analyses, and information theoretical approaches have been used for network inference from population-level data. Combining multiple inferred networks to form a community or ensemble prediction often confers slight but consistent improvements in the quality of the predicted network, but how to best combine and weight different methods to form a consensus prediction is poorly understood. Given our current understanding, a reasonable approach to generate such ensemble predictions would be to include information derived using different classes of inference algorithms (since these are known to show different biases; ) but also to ensure that within each class, we develop the best-performing algorithm based on a given statistical methodology.

4More recently, the increasing availability of single-cell expression data has led to the development of several computational and statistical approaches aimed at gaining new insight into cell fate decisions and transitions between cell states. Identifying the associated changes in transcriptional state and regulatory interactions that contribute to controlling these processes are key aims of many single-cell transcriptomic studies. Several pseudotemporal ordering algorithms have been developed that aim to place cells in an inferred temporal order based on similarities in their transcriptional states, since true single-cell temporal data (where large numbers of genes are assayed) are not feasible to collect at present. These methods often make strong assumptions about developmental processes (which have been questioned; ) and, in most cases, uncertainties in the inferred order are likely to affect and bias downstream analyses. Network inference methods, in contrast, explore statistical dependencies between genes (from the observed distributions of expression levels across a given population of cells) and identify those that may be indicative of functional relationships, without necessarily making such strong assumptions about the nature of cell transitions; each inferred edge is a hypotheses that can be tested. Information theoretical measures, in particular, are particularly parsimonious in the assumptions that they make compared with, e.g., ordinary differential equation (ODE)-based regression approaches or simple graphical models. However, network inference using single-cell data (Figure 1) remains relatively unexplored with only a few notable examples (e.g.,,, ) extending beyond simple, potentially simplistic, correlation-based analyses, perhaps due to the difficulties of directly applying and interpreting the results of existing methods designed to deal with population-level data.

5Single-cell data are notoriously complex and present new challenges for statistical analysis: technical noise is difficult to distinguish from genuine biological variability, the relative contributions and impact of different sources of noise are poorly understood, and numerous factors contribute to the biological heterogeneity observed within cell populations. However, these data also offer substantial advantages over population-level data that we can exploit in order to learn about the structure of GRNs governing the observed changes in gene expression. Firstly, datasets are large, routinely comprising expression measurements from hundreds or thousands of cells, and continuing advances in single-cell technologies will allow further increases in sample sizes. In addition, single-cell data inherently provide the variability required to detect statistical relationships between gene expression profiles (interpreted as putative functional relationships), whereas population-based studies need to introduce this variability by observing cell populations under different conditions, e.g., using time-series or perturbation studies. While we can collect single-cell time-series data, even data collected at one time point will contain variability due to (1) asynchrony of cells within a population (in terms of progression through a biological process), and (2) biological heterogeneity and often the presence of multiple cell (sub)types.

6Here, we introduce an information-theory-based GRN inference algorithm designed to make use of these features of single-cell data. Information theory provides a set of measures, chiefly among them mutual information (MI), that allow us to characterize statistical dependencies between pairs or groups of random variables without making assumptions about the nature of the dependencies. MI has considerable advantages over simpler measures such as (Pearson) correlation, as it is capable of capturing complex non-linear and non-monotonic dependencies, and reflecting the dynamics between pairs or groups of genes. Calculating MI involves estimating pairwise joint probability distributions, generally requiring density estimation or data discretization, and the accuracy of these estimates depends on the sample sizes. Single-cell datasets are sufficiently large to allow us to accurately estimate probability distributions between more than two variables, and thus capitalize on recent developments in multivariate information (MVI) theory. Quantifying the information between three or more variables has long been problematic, and the most widely used measures, such as interaction information, are known to have serious flaws. The recently introduced partial information decomposition (PID) both explains and solves its predecessors' problems and provides a meaningful measure of MVI; PID and other measures are described in detail in Box 1. Our algorithm uses PID to analyze the statistical relationships between triplets of variables to generate undirected networks highlighting putative functional interactions between genes.The entropy, H(X), quantifies the uncertainty in the probability distribution, p(x), of a random variable X. For a discrete random variable, the entropy is given bywhich is maximal for a uniform distribution. If we consider the mRNA expression level of a gene to be the variable X, then a gene that is expressed differently across a set of cells will have a higher entropy than a more consistently expressed gene. When considering the relationship between X and a second random variable, Y, we quantify the information that one variable provides about the other using the MI:This quantifies the difference between the joint entropy, H(X,Y), and the joint entropy assuming independence of X and Y, H(X) + H(Y), and thus provides a non-negative, symmetric measure of the statistical dependency between the two variables. For a pair of genes with more co-ocoordinated expression, their observed joint entropy will be lower, and hence they will have higher mutual information. Given a third variable, Z, the conditional mutual information (CMI),quantifies the information between X and Y given knowledge of Z. This tells us the extent to which knowing the expression of one gene additionally informs us about the expression of a second gene, given that we already know the expression of a third gene.A number of MVI measures have been defined that aim to quantify the statistical dependencies between three or more variables, but there is little consensus as to the most appropriate metric. Arguably, the most widely used is interaction information, which for three genes quantifies the extra information between any two of the genes X, Y, and Z, when the third is known compared with when it is not known:Interaction information has received much criticism since it can be (1) zero between dependent variables, when MI and CMI are equal but non-zero, and (2) negative, when MI is greater than CMI. The problem is that, despite its name, interaction information is not a quantity of information between a set of variables; rather it quantifies the balance between MI and CMI. In fact, MVI is difficult to summarize with a single quantity, because there are different ways in which information can be shared by three or more variables.An important recent development in information theory is the introduction of PID. In the three-variable case, PID considers the information provided by a set of source variables (or genes), S = {X, Y}, about another target variable, Z, partitioned into redundant, synergistic, and unique information. Redundant information is the portion of information about Z that can be provided by either variable in S alone; the unique contribution from X (or Y) is the portion of information provided only by X (or only Y); and the synergistic information is the portion of information that is only provided by knowledge of both X and Y. Thus, the PID between the set S and the target variable Z is equal to the sum of the four partial information terms,where UniqueY(Z;X) is the unique information between source variable X and target variable Z when the other source variable is Y. In the context of gene expression, the terms source and target do not imply any mechanistic assumptions; rather, they refer to quantifying to what extent knowledge of the source genes informs us about the target gene.To calculate the PID terms, the redundant information is first calculated using the specific information, Ispec, which quantifies the information provided by one variable about a specific state of another variable. The “state” of a gene in a given cell refers to which discrete bin its mRNA level falls within, once we have discretized the expression data. If we consider the information provided by X about state z of variable Z,If we consider S = {X, Y} and a target variable Z, the redundant contribution is calculated by comparing the amount of information provided by each variable within S about each state of the target Z,The unique information terms can be calculated from the redundant information and the pairwise MI, via the relationship,In other words, the pairwise MI between two variables can be partitioned into a redundant and unique component, given a third variable; it is this relationship that we exploit further in the Results section. We note that although PID is not a symmetric measure, UniqueY(Z;X) + Redundancy (Z;X, Y) will be equal to UniqueY(X;Z) + Redundancy(X;Z,Y) because MI is symmetric.Finally the synergistic information can be calculated via the interaction information (Equation 4), which turns out to be the difference between the synergistic and redundant information:

7The entropy, H(X), quantifies the uncertainty in the probability distribution, p(x), of a random variable X. For a discrete random variable, the entropy is given bywhich is maximal for a uniform distribution. If we consider the mRNA expression level of a gene to be the variable X, then a gene that is expressed differently across a set of cells will have a higher entropy than a more consistently expressed gene. When considering the relationship between X and a second random variable, Y, we quantify the information that one variable provides about the other using the MI:

8This quantifies the difference between the joint entropy, H(X,Y), and the joint entropy assuming independence of X and Y, H(X) + H(Y), and thus provides a non-negative, symmetric measure of the statistical dependency between the two variables. For a pair of genes with more co-ocoordinated expression, their observed joint entropy will be lower, and hence they will have higher mutual information. Given a third variable, Z, the conditional mutual information (CMI),quantifies the information between X and Y given knowledge of Z. This tells us the extent to which knowing the expression of one gene additionally informs us about the expression of a second gene, given that we already know the expression of a third gene.

9A number of MVI measures have been defined that aim to quantify the statistical dependencies between three or more variables, but there is little consensus as to the most appropriate metric. Arguably, the most widely used is interaction information, which for three genes quantifies the extra information between any two of the genes X, Y, and Z, when the third is known compared with when it is not known:

10Interaction information has received much criticism since it can be (1) zero between dependent variables, when MI and CMI are equal but non-zero, and (2) negative, when MI is greater than CMI. The problem is that, despite its name, interaction information is not a quantity of information between a set of variables; rather it quantifies the balance between MI and CMI. In fact, MVI is difficult to summarize with a single quantity, because there are different ways in which information can be shared by three or more variables.

11An important recent development in information theory is the introduction of PID. In the three-variable case, PID considers the information provided by a set of source variables (or genes), S = {X, Y}, about another target variable, Z, partitioned into redundant, synergistic, and unique information. Redundant information is the portion of information about Z that can be provided by either variable in S alone; the unique contribution from X (or Y) is the portion of information provided only by X (or only Y); and the synergistic information is the portion of information that is only provided by knowledge of both X and Y. Thus, the PID between the set S and the target variable Z is equal to the sum of the four partial information terms,where UniqueY(Z;X) is the unique information between source variable X and target variable Z when the other source variable is Y. In the context of gene expression, the terms source and target do not imply any mechanistic assumptions; rather, they refer to quantifying to what extent knowledge of the source genes informs us about the target gene.

12To calculate the PID terms, the redundant information is first calculated using the specific information, Ispec, which quantifies the information provided by one variable about a specific state of another variable. The “state” of a gene in a given cell refers to which discrete bin its mRNA level falls within, once we have discretized the expression data. If we consider the information provided by X about state z of variable Z,

13If we consider S = {X, Y} and a target variable Z, the redundant contribution is calculated by comparing the amount of information provided by each variable within S about each state of the target Z,

14The unique information terms can be calculated from the redundant information and the pairwise MI, via the relationship,

15In other words, the pairwise MI between two variables can be partitioned into a redundant and unique component, given a third variable; it is this relationship that we exploit further in the Results section. We note that although PID is not a symmetric measure, UniqueY(Z;X) + Redundancy (Z;X, Y) will be equal to UniqueY(X;Z) + Redundancy(X;Z,Y) because MI is symmetric.

16Finally the synergistic information can be calculated via the interaction information (Equation 4), which turns out to be the difference between the synergistic and redundant information:

17In this paper, we describe an inference algorithm based on the MVI measure PID and use extensive in silico analyses to demonstrate (1) the consistent improvement over existing algorithms based on pairwise MI, and (2) the suitability of this method for analyzing single-cell data, before illustrating its application to several experimental datasets. Such in silico analyses are critical for quantitatively assessing network inference approaches as, unlike with real biological systems, we have knowledge of the “true” GRN underlying the observed data. These results demonstrate that the larger sample sizes of single-cell data are vital for our method and that they profoundly improve the performance of information-based methods in general. We thoroughly explore the factors that influence algorithm performance — in particular the choice of discretization algorithms and probability distribution estimators — in order to provide evidence-based guidelines for the use of information-theory-based methods for network inference. We emphasize the importance of considering the different sources of heterogeneity within single-cell data so that we can take advantage of the variation of interest — e.g., that associated with progression through a biological process such as differentiation. Our examples using experimental data demonstrate how our inference method can be combined with existing computational and statistical methods (e.g., clustering and dimensionality reduction) to infer networks from carefully chosen subsets of single-cell data in order to address particular questions about cellular processes. We consider three single-cell transcriptomic datasets here and additionally refer to a related manuscript ( [this issue of Cell Systems]), in which we use our framework to infer changing regulatory subnetworks over the course of neural progenitor development in mouse embryos, and to suggest candidate genes for maintaining cellular states and driving state transitions. Finally, we provide a fast, open-source implementation of our methods to enable easy application to other single-cell datasets.

18We first investigate the usefulness of PID for inferring network edges using data generated from in silico models. We use stochastic simulations from simple directed 3-node networks of varying topologies and estimate PID values (redundant, synergistic, and unique information, defined in Box 1) from these simulated data. Simulations were generated from two model definitions, based on thermodynamic or mass-action kinetics, both commonly used in systems biology to represent gene regulation, as described in STAR Methods. A distinctive pattern is apparent in networks with a single directed edge between two genes (“one-edge” topology; Figure 2A) — the unique information between the two connected genes is notably higher than both the unique information between unconnected genes and the redundancy values between all three genes. With increasing numbers of edges within the network, this pattern is lost; this makes sense intuitively as with higher connectivity we expect to see increased synergistic or redundant contributions (Figures 2A and S1). Note that the pattern can only be observed under simulation conditions that generate variability in the observed variables (i.e., statistical relationships are not detectable when the system is at steady state; Figure S1).

19To explore whether this pattern also occurs for triplets of nodes embedded in large networks, we consider time-series expression data simulated from five different 50-gene networks generated by GeneNetWeaver. This software generates stochastic simulations from dynamical models that represent transcription and translation using a thermodynamic approach, with network structures that are inspired by known gene connectivity patterns in Escherichia coli and Saccharomyces cerevisiae (Figure 2B), and it has become a standard tool for performance evaluation of network inference algorithms. PID values are estimated from these data for every triplet of genes within the networks, and each triplet is classified according to its topology — six topological arrangements are possible given the model assumptions (maximum of one edge between each pair of nodes, no self-regulation, and no feedback loops). Mean PID values are calculated for each group, and the same distinctive pattern (high unique versus redundant contributions for connected genes) is apparent for triplets with a single directed edge (Figure 2C). As with the 3-node simulations, the pattern is lost in topologies with more connections; however, the “unconnected” and “one-edge” topologies are by far the most prevalent in all networks — comprising 64.4%–93.2% and 6.3%–29.7% of all triplets, respectively, and jointly comprising over 90% of triplets in every network (Table S1 and Figure S2).

20Examining in silico data from a 10-gene network suggests that the relative size of the unique information compared with the redundancy — i.e., the proportion of MI accounted for by the unique contribution (Equation 8 in Box 1) — is more informative than the absolute unique information (Figure 3A). Using data from the 50-node networks, we confirm that this is the case more generally. We take all triplets of the unconnected and one-edge topologies (with each gene being able to take part in multiple triplets); then for each triplet, treating each gene in turn as the target, we estimate (1) the redundant information between all three genes, (2) the unique information between one source gene and the target, and (3) the unique information between the second source gene and the target. We plot the ratio of unique information to MI (the sum of unique and redundant information) for each pair of genes in each triplet and find that this ratio is higher in general between connected pairs (Figure 3B).

21In a network of n genes, given a pair of genes X and Y, there are n − 2 gene triplets involving the pair. The MI between X and Y, I(X;Y), is unaffected by the choice of the third gene, Z, because MI is a pairwise measure, but the unique information between X and Y, UniqueZ(X;Y), varies depending on Z. Furthermore, the difference between I(X;Y) and UniqueZ(X;Y) is equal to the redundancy between all three genes (Equation 8), meaning that we can regard the ratio UniqueZ(X;Y)/I(X;Y) as capturing the proportion of MI that is accounted for by unique information between X and Y, as opposed to redundant information between all three genes. We note that the vast majority of all possible gene triplets in our E. coli and S. cerevisiae networks are of the unconnected or one-edge topology (Table S1), and that the ratio UniqueZ(X;Y):I(X;Y) is higher between connected pairs in a one-edge triple (Figure 3B). Therefore, we would expect that if X and Y are connected, then most of the triplets made with X, Y, and all n − 2 possible Z in turn, would be of the one-edge topology — and likewise if X and Y are unconnected, then most triples involving X and Y would be of the unconnected topology — so UniqueZ(X;Y)/I(X;Y) would in general be higher if X and Y were connected.

22We define the proportional unique contribution (PUC) between two genes X and Y as the sum of this ratio calculated using every other gene Z in a network (where S is the complete set of genes):

23this measure may be thought of as capturing the mean proportion of MI between two genes X and Y that is accounted for by the unique information. Note that the PID unique measure is not symmetric, so for each pair of genes we consider each as the target in turn (hence we include both UniqueZ(X;Y) and UniqueZ(Y;X) terms in Equation 10).

24In our network inference algorithm (Figure 4A), the redundancy and unique information contributions are first estimated for every gene triplet, then the PUC is calculated for each pair of genes in the network (Equation 10). Finding a threshold for defining an edge at this stage is problematic, because the distributions of PUC scores differ between genes (see Figure 4B), thus setting a global threshold for PUC scores across the whole network risks biasing the results by factors such as expression variability. This was previously observed with MI and led to the development of measures that take into account the network context, central to the context likelihood of relatedness (CLR) algorithm. A similar solution is employed here: an empirical probability distribution is estimated from the PUC scores for each gene, and the confidence of an edge between a pair of genes is given bywhere FX(U) is the cumulative distribution function of all the PUC scores involving gene X (here, we assume either a Gamma or Gaussian empirical probability distribution). This effectively identifies the most important interactions per gene, rather than just taking the highest pairwise scores across the whole network.

25We compare our algorithm against several common information-theory-based network inference methods, and thus briefly summarize these existing approaches in Box 2. All methods start from the pairwise MI matrix and then use it in different ways. Even compiling the MI matrix is, however, fraught with potential problems: the manner in which data are treated (e.g., discretization) and the estimator used for the entropy and MI both affect the performance of the algorithms. When comparing different approaches, it is therefore important to ensure that discretization and estimation of MI are done identically. Without this, it becomes impossible to disentangle the relative strengths and weaknesses of the different approaches that are based on and interpret MI values. In the STAR Methods, we discuss the different estimators and discretization approaches that we use and which are implemented in the InformationMeasures.jl package (see Software section). Our comparisons with existing methods thus always start from the same MI matrix.Relevance networks use the MI estimates (or, in some cases, correlation) in order to detect edges. As there is no reliable universal way of determining the statistical significance of MI values, a threshold is typically chosen to determine which edges are present. This fails to account for the fact that MI may be increased for nodes X and Z even though they only indirectly interact via an intermediate node Y. The Data Processing Inequality (DPI) allows us to sort out some of these cases by virtue of the relationshipwhich holds whenever X, Y, and Z form a Markov Chain. Post-processing of the MI values using the DPI is at the core of the popular ARACNE algorithm. Thresholds on the pairwise MIs are used to identify likely dependent pairs X, Z; MI values above the threshold are then considered with every possible other node Y in light of the DPI.Given that MI values are affected by a number of factors, including especially the variability of each individual random variable, any global a priori threshold may be highly problematic: it will give rise to false positives as well as false negatives. In the CLR algorithm, the MI between X and Z is considered against all MI values for pairings of X and Z with all other variables Y. Thus, the threshold for each pair will reflect the variabilities of both genes, as well as their relative levels of statistical dependence on other genes. MRNET aims to identify a minimally redundant but maximally explanatory set of variables/predictors for each target gene X in a greedy manner.There have been attempts at using conditional mutual information (Equation 3), interaction information (Equation 4), or related concepts, for network inference (although not applied to single-cell data). These would have to deal with the known difficulties of interpretation (which do not arise in relation to PID-based measures), which may explain the lack of their widespread uptake.

26Relevance networks use the MI estimates (or, in some cases, correlation) in order to detect edges. As there is no reliable universal way of determining the statistical significance of MI values, a threshold is typically chosen to determine which edges are present. This fails to account for the fact that MI may be increased for nodes X and Z even though they only indirectly interact via an intermediate node Y. The Data Processing Inequality (DPI) allows us to sort out some of these cases by virtue of the relationshipwhich holds whenever X, Y, and Z form a Markov Chain. Post-processing of the MI values using the DPI is at the core of the popular ARACNE algorithm. Thresholds on the pairwise MIs are used to identify likely dependent pairs X, Z; MI values above the threshold are then considered with every possible other node Y in light of the DPI.

27Given that MI values are affected by a number of factors, including especially the variability of each individual random variable, any global a priori threshold may be highly problematic: it will give rise to false positives as well as false negatives. In the CLR algorithm, the MI between X and Z is considered against all MI values for pairings of X and Z with all other variables Y. Thus, the threshold for each pair will reflect the variabilities of both genes, as well as their relative levels of statistical dependence on other genes. MRNET aims to identify a minimally redundant but maximally explanatory set of variables/predictors for each target gene X in a greedy manner.

28There have been attempts at using conditional mutual information (Equation 3), interaction information (Equation 4), or related concepts, for network inference (although not applied to single-cell data). These would have to deal with the known difficulties of interpretation (which do not arise in relation to PID-based measures), which may explain the lack of their widespread uptake.

29Undirected networks are inferred from in silico datasets (described in the STAR Methods) for five 50-gene networks and five 100-gene networks using ARACNE, CLR, MI (relevance networks), MRNET, and the PID-based algorithm, PID and context (PIDC). We also include results for the raw PUC score, without the network context step. Accuracy of the inferred networks is evaluated using the area under the precision-recall curve (AUPR), rather than the receiver operating characteristic curve (AUROC), which is inappropriate for judging network inference methods as real networks are typically sparse; see STAR Methods for definitions and a more detailed discussion.

30PIDC performs favorably compared with the other algorithms (Figure 5), particularly in the larger networks. The raw PUC score outperforms the raw MI score, indicating the value of higher-order information; and CLR outperforms the other MI-based approaches, indicating the value of network context, in agreement with the previous comparisons (Figures 5A and 5C). This effect is robust to simulated technical noise (Figure 5B) and becomes more evident the larger the number of “cells” in the dataset (Figure 5A); as real data become more accurate and sample sizes increase, we expect the performance of PIDC/PUC to improve as the estimation of three-dimensional dependencies becomes more accurate. Also, unlike CLR, it was designed to capture higher-order dependencies and to distinguish between direct and indirect interactions.

31The extensive analyses using simulated data are necessary to validate our algorithm and provide quantitative comparisons with existing methods. When working with real experimental datasets, we of course do not know the true underlying network, and thus rely on identifying relationships that are consistent with our current biological knowledge about the systems we are studying. Here, we apply our algorithm to three published experimental datasets and, in a related manuscript, illustrate how it can be used as part of a thorough modeling analysis of neuronal differentiation of mouse embryonic stem cells.

32used single-cell qPCR (sc-qPCR) to study megakaryocyte-erythroid progenitors (MEP) during human hematopoiesis. Their analysis revealed the existence of subpopulation structure in this class of cells: two groups of cells are primed preferentially for a particular cell fate — megakaryocytic (MK-MEP) or erythroid (E-MEP) — while a third group of multi-potent progenitors (Pre-MEP) retain some myeloid differentiation capacity (Figure 6A). Here, we apply our PIDC algorithm to their complete dataset and infer a candidate network that depicts statistical dependencies among the analyzed genes (Figure 6B). Given we are interested in genes involved in differentiation processes, we also consider networks inferred using overlapping subsets of the data and color each edge in Figure 6B according to their presence in these additional networks. Edges that are present in the network constructed using Pre-MEP and E-MEP cells, but not that based on Pre-MEP and MK-MEP cells, are colored red (i.e., erythroid-specific edges), and the reverse (megakaryocytic specific) are shown in blue (edges present in both, or just in the original complete network are shown in gray). Consistent with existing knowledge about these two lineages, we see a cluster dominated by blue edges that comprises known megakaryocytic genes (e.g., CD9, LOX, VWF, NFIB, CD61, TGFB1) and a cluster with several red edges comprising known erythroid genes (e.g., CD36, KLF1, LEF1, CNRIP1, TMOD1, ANK1, DHRS3). Networks based on pairwise MI scores from the same data show skewed degree distributions with many nodes unconnected (Figures 6C and 6D).

33We next consider an sc-qPCR dataset comprising expression measurements of selected genes during early embryonic development (from oocyte to 64-cell blastocyst stages); Figure 6E shows the resulting inferred network. We again infer additional networks using subsets of these data — here we use overlapping subsets of “early” and “late” cells to reveal any temporal dependencies in the detected interactions. A number of known relationships between genes are apparent in the network, e.g., upregulation of Cdx2 and Gata3 transcription factors (TFs) during the 8-cell to morula transition is identified as an edge in the early network; while the co-expression of primitive endoderm specific TFs Creb312 and Sox17 is detected as an edge in the late network (consistent with the appearance of distinct cell types, including primitive endoderm cells, in the blastocyst). A cluster of known pluripotency and reprogramming factors is also identified in the network (Pou5f1, Nanog, Esrrb, Klf2, and Klf4) — Sox2, another key reprogramming factor, is not connected with these genes but is known to be up-regulated later than the other factors.

34studied embryonic hematopoietic development and used sc-qPCR data to develop a Boolean network model of the GRN underlying blood development. In Figure 6F, we compare the networks inferred from these same data using PIDC and MI — we find the inferred PIDC network shares a higher number of edges with the Boolean network model than the network constructed using MI values alone. Although we, of course, do not know the true structure of the GRN in this case (this is only feasible when using in silico data), the Boolean model was shown to capture key cell states observed experimentally and generated several experimentally validated predictions, thus we use this as a benchmark to indicate the biological plausibility of our inferred networks.

35In a related companion paper, we apply our PIDC algorithm to sc-qPCR data collected from cells undergoing differentiation from a pluripotent ground state toward a committed neuronal lineage, via a primed epiblast-like state. A total of 547 cells were sampled at seven time points spanning 7 days; we analyze expression measurements of 74 genes, including known regulators of pluripotency and neuronal differentiation (see for details). We first assign the cells to three robust groups that correspond with developmental stage (using k-means clustering) and then infer networks using data from all the cells or from overlapping subsets containing cells at earlier or later stages of development. Comparing the networks obtained using different subsets of cells allows us to observe any temporal dependencies in the inferred interactions. Using an unsupervised community detection algorithm, we find that the network of inferred (co-)regulatory relationships contains several communities (or modules) of genes displaying high connectivity within each community. Three of these communities show distinct temporal dependencies in connectivity and comprise genes known to play roles at different stages of differentiation. Our analysis thus identifies modules of genes that undergo coordinated changes in expression as cells progress through development and putative gene interactions that may be involved in regulating these transitions in cell state.

36Any comparative analysis of information-based GRN inference algorithms is influenced by a number of decisions, in particular: (1) how the data are discretized, (2) the choice of MI estimator, and (3) the metric used to evaluate performance. We discuss the impact of each of these decisions and offer guidelines for future analyses, before discussing the use of single-cell RNA sequencing (scRNA-seq) datasets.

37The information measures described here all rely on estimates of discrete probability distributions. Normalized mRNA expression data are generally continuous, but estimating the distributions for continuous random variables is fraught with problems. Several algorithms and heuristics have been developed to discretize data and estimate empirical probability distributions from the resulting discrete frequencies. We investigated two methods for discretization, along with four MI estimators, as described in STAR Methods.

38All estimators produce fairly accurate estimates of joint entropies for up to two uniformly distributed random variables, but in higher dimensions, performance varies according to the distribution (Tables S2 and S3), making it difficult to identify the most appropriate estimator for experimental data. Rank agreement between the estimators is good when the data are discretized using Bayesian blocks, however, diminishing the importance of the choice of estimator (Figure S4). In light of these findings we advise using Bayesian blocks (an adaptive discretization algorithm that allows variable-width bins), and since the underlying distribution is usually unknown, we favor the maximum likelihood estimator due to its simplicity.

39The choice of discretization method and estimator influences the performance of the inference algorithms, with effects varying depending on the algorithm and on the true network (Figure S3). Sampling frequency and dataset size also have an effect (Figure 5A), with performance increasing in line with dataset size and decreasing in larger networks and with a higher rate of technical dropout errors (Figure 5B). Due to the number of influential factors, we advise caution when interpreting the results of this or any such comparison as an exhaustive exploration of these factors is not feasible; however, we note that the PIDC algorithm performs well in general across the many tested combinations of discretization methods, estimators, and datasets.

40The metric used to evaluate algorithms also affects their apparent performance, evident here in the higher scores for AUROC than AUPR (Figure S3B). This is a well-documented phenomenon, caused by the true negatives (unconnected node pairs) in a GRN vastly outnumbering the true positives (edges); for example, the E. coli 1 100-gene network contains 125 edges and 4,825 unconnected pairs (Table S1). AUROC equally rewards the prediction of an edge and a non-edge, meaning that the score for any algorithm that mostly (or even exclusively) predicts non-edges will be misleadingly inflated, however well or badly it predicts edges. AUPR is therefore the more meaningful measure, despite AUROC being widely used.

41Here, we have illustrated the application of our method using several single-cell qPCR datasets; scRNA-seq experiments generate much larger datasets comprising expression measurements for thousands of genes. When analyzing these data, a subset of (up to hundreds of) genes should first be selected — both to make the network inference analysis computationally tractable (see Software section) but also to aid in interpreting the results. There are many potential approaches to selecting gene subsets, depending on the purpose of the analysis and the level of existing knowledge about the specific system being studied. We may wish to select genes likely to be involved in the process of interest based on prior knowledge (and/or functional annotations of genes) such as known TFs, similar to the way that genes are selected for analysis in qPCR experiments. However, we can also make use of the huge array of statistical and computational methods that have been developed to analyze scRNA-seq experiments and select subsets of genes based on their observed gene expression patterns. For example, we could select those showing higher than expected levels of expression variability, those that show differential expression between cell states or over time (using the results of clustering or pseudotemporal ordering algorithms), or cluster genes by similar expression profiles and select representative genes from each cluster. As is the case for other network inference algorithms, genes with no variability in mRNA expression are uninformative and should always be removed prior to analysis. In addition, due to the prevalence of zeros and lack of sensitivity of single-cell experiments, many genes (particularly those expressed at low levels) will not be reliably detected so we can also exclude those without detectable expression in a large proportion of cells.

42A new open-source package for estimating MVI measures is implemented in the Julia programming language. The package, named InformationMeasures.jl supports information measures such as entropy, MI, conditional mutual information (CMI), and PID; the maximum likelihood, Miller-Madow, Dirichlet, and shrinkage estimators; and the Bayesian blocks, uniform width, and uniform count discretization methods.

43Julia was chosen for its speed (Figure 7A), clear mathematical syntax, growing availability of libraries, and good integration with other languages. The existing Discretizers.jl package is used to implement the discretization methods; in some of our initial analyses we used the AstroML Python implementation of the Bayesian blocks algorithm. In order to meet a wide range of requirements, the package can be used simply for discretizing data or to calculate information measures using pre-discretized data or probability distributions that have been estimated elsewhere.

44A Julia package implementing PIDC and other inference algorithms is available, along with tutorials and our simulated datasets at https://github.com/Tchanders/network_inference_tutorials. Our algorithm has complexity O(n3) in the number of genes, but the speed of our Julia implementation means that inference time is comparable with widely used implementations of the lower-complexity algorithms (Figure 7A). The complexity in the number of cells depends on the discretization method: the recommended Bayesian blocks method scales less well than the uniform width method but nevertheless produces results for several thousands of cells for a network of hundreds of genes within a practically useful timescale (Figures 7B and 7C).

45Here, we have introduced a network inference algorithm based on PID, an easily interpretable MVI measure that allows us to explore statistical dependencies between multiple genes in detail. Our PIDC algorithm identifies putative functional relationships between genes based on the unique contribution to pairwise MI (Equation 8) combined with information about the local network context of each gene. We use extensive performance comparisons (Figure 5) to demonstrate the value of using both higher-order information measures and network context and to illustrate that the large sample sizes provided by single-cell data are critical to our algorithm's success. Like other studies, we find that the methods chosen to discretize data and estimate entropies and probability distributions affect algorithm performance considerably (Figure S3) — too often, the impact of these choices has been ignored. A fast, open-source software package provides an easy way for users to explore such factors when applying our method.

46Although single-cell data have many potential advantages over bulk transcriptomic data for network reconstruction approaches — particularly sample size, inherent variability, and ability to detect subpopulation structure — they also pose notable challenges. Technical noise and other biological sources of heterogeneity (e.g., transcription stochasticity, other cellular processes) can, to some extent, impede our ability to detect informative statistical dependencies; there are clear theoretical benefits of using information theory methods when we expect the observed dependencies to be complex and non-linear. The relative contributions of different sources of noise and variation are still poorly understood, and good noise models are lacking. We therefore relied on a simple model of dropout events to mimic the zero-rich nature of these data in our simulations — as expected, technical noise reduces our ability to recover the true network structure (Figure 5B). The best-performing algorithms, PIDC and CLR, both aim to capture the most important interactions for each node in turn, rather than the highest dependencies across the whole dataset. This leads to well-connected inferred graphs (e.g., Figure 6C) but can also help to address the influence of potential confounding factors (such as the cell cycle) when working with single-cell data. Variation in cell-cycle stage causes large-scale changes in cell transcriptional states — and we would expect it to induce stronger statistical dependencies among the affected genes. However, using empirical distributions to take network context into account (as in PIDC and CLR) will at least partly mitigate the influence of any such confounding factors.

47Integrating our method with other single-cell analyses allows us to select subsets of cells and genes that are most informative about our specific biological questions. Firstly, there are many sources of biological heterogeneity in single-cell data so it is important to focus on the variation of interest, e.g., when studying developmental processes, we should aim to analyze collections of data where we expect cell differentiation to be the major source of variation. We can use established methods for analyzing the subpopulation structure of single-cell data — such as clustering, dimensionality reduction, and pseudotemporal ordering algorithms — e.g., to identify distinct cell subtypes or alternative differentiation pathways. Here, we have shown how this allows us to focus on functional relationships involved in developmental transitions and the relative timing of transcriptional changes. For example, in, we use clusters of developmentally similar cells to examine how the activity of (co-)regulatory relationships changes during development, making initial steps toward defining mechanistic models of neural differentiation — the network structure and temporal information allow us to suggest candidate genes for maintaining cell states or driving developmental transitions. Although scRNA-seq generates data for thousands of genes, we recommend selecting meaningful subsets for network inference — either based on biological knowledge or by using gene expression variability and patterns. Data imputation methods (e.g., ) or average expression measurements over small groups of similar cells may help address the challenges of noise and low coverage.

48For any statistical approach, it is important to consider the potential limitations. Many are general to network inference approaches aiming to reconstruct GRN structure from mRNA profiles and are discussed in depth elsewhere, but it is worth emphasizing a few key points that affect how we should interpret and use our results. Firstly, we can only detect relationships where there is sufficient variability in gene expression observed under the chosen experimental conditions. Functional interactions are only detectable if they induce changes in transcriptional state that persist over a reasonable timescale — we will not, for example, detect rapidly fluctuating changes, as the transient changes in mRNA levels will not result in observed statistical dependencies across cells. As well as functional regulatory relationships, we are likely to also identify co-regulatory relationships where genes under the influence of the same regulator show coordinated expression changes. Without making further assumptions, or using perturbation or temporal data, we cannot distinguish causal relationships; however, in many instances it will still be informative to learn which sets of genes respond in a coordinated manner. In cases where the assumptions made by pseudotemporal ordering algorithms are justified, we can potentially use this information to infer causality and directionality of gene interactions. It is of course unrealistic to expect every cell to follow precisely the same route through transcriptional space; we do instead make the modest assumption that there are key changes in transcriptional state that must occur in order for cells to respond appropriately to environmental and developmental cues, and that these will be subject to conserved regulatory mechanisms.

49Methods for exploring high-throughput single-cell datasets and identifying putative functional relationships between genes are clearly needed. As with all network inference methods, we cannot expect to reconstruct the exact structure of the underlying biological networks but instead view such methods as tools to explore the data; generate hypotheses; represent the current state of understanding; and guide further experiments, model development, and analyses. Validating or invalidating these hypotheses experimentally may, of course, lead to revised network models — like any mathematical model, it should be subject to refinement as new facts are assembled and new insights are gained.

50Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Ann C. Babtie (a.babtie@imperial.ac.uk).

51In order to use the entropy estimators described here, continuous datasets must first be discretized. A number of algorithms exist to define the total number and boundaries of the resulting partitions (or bins). One common simple approach is to use bins of equal width, with the number of bins determined heuristically; e.g., here, we use the nearest integer to the square root of the size of the dataset,. A more sophisticated approach is the Bayesian blocks algorithm, in which the number and widths of bins are chosen by optimizing a fitness function, without constraining the bins to be of equal width.

52When a dataset is large enough, the empirical frequencies can be considered to be an approximation of the true probabilities, referred to here as the maximum likelihood approach. For sparser datasets, a number of methods have been developed either for estimating the probability distribution from a set of frequencies, such as the Dirichlet estimator and the shrinkage estimator, or for estimating the entropy directly, as in the Miller-Madow estimator.

53The Dirichlet estimator refers to a group of Bayesian estimators that take a Dirichlet distribution as prior but each with different parameters. There is no consensus on the best parameters to use, despite several proposed alternatives; here, we use the same parameter, 1, for each bin unless otherwise stated.

54The shrinkage estimator is also Bayesian, compromising between the observed frequencies, unbiased but with a high variance, and a prior (or target) distribution, biased but with low variance. The estimate is affected by both the choice of target distribution and the weight given to the target (or shrinkage intensity). In the current analysis, the optimal shrinkage intensity is calculated as described in, and the target distribution is the uniform distribution.

55The Miller-Madow estimator is an entropy bias correction that does not estimate the probability distribution, and therefore cannot be meaningfully applied to higher-order information measures. Despite this, it has been applied for the comparison of different MI-based algorithms, and so it is included in this analysis, with the caveat that its meaning is unclear.

56We considered six three-gene topologies (Figure S1) and used the Gillespie algorithm to generate stochastic simulations of gene expression time-course data using two alternative model definitions (based on thermodynamic or mass-action kinetics). In both cases, we include an additional activating stimulus that is present from halfway through the simulation time course. This additional stimulus acts to perturb the system away from a steady state, driving changes in gene expression that are necessary for relationships between genes to be observable. This section describes these model and simulation details.

57The first thermodynamic model includes seven species: mRNA (xi) and protein (yi) corresponding to three genes (i = 1, 2, 3), and a stimulating ligand (s) which targets a selected gene. We define the following reaction types with associated propensities:to represent mRNA decay, protein decay, translation, and transcription, respectively, where λ is the protein/mRNA decay rate, translation is modeled according to saturation kinetics (with maximum rate αtranslation), i.e.,

58and transcription rates depend on the concentration of any regulating proteins (including the stimulating ligand s if present) according to the relationshipwhere αtranscription is a constant transcription rate, M is the total number of possible states Sm for gene i (either unbound or bound by one or two regulating proteins), and αm is the relative activation rate for each state. The probability of each state, P{Sm}, depends on the concentrations of the regulating proteins, modeled according to standard thermodynamic principles (see e.g., for details). For example, if a gene has two possible regulators (proteins yj,yk), we calculate the mean activation of transcription of the target gene i using the function,where χj = (yj/kj), kj is the dissociation constant, and the possible states of the gene are unbound, bound by yj or yk alone, or by both yj and yk. For our models, the maximum number of regulators for a given gene is three (proteins yj and yk plus the stimulating ligand s), but we assume a maximum of two regulators can bind the gene at any one time (so we consider the gene states with each possible pair of ligands bound but not a state with all three bound).

59For the thermodynamic model, we simulated time courses from time 0 to 1,000, with the stimulating ligand s present from time 500 at a constant level of 20 molecules, and recorded the system state at 41 equally spaced intervals. We repeated the simulations 25 times — and calculated PID scores using the resulting data (i.e., 1,025 data points or cells were used to calculate each PID measure)— with the stimulus targeting each of the three genes in the network in turn (75 simulations in total, 3 sets of PID scores). We randomly sample initial mRNA (xi) and protein (yi) levels from a U(0,5) distribution; we perform simulations for five different initial conditions and plot the mean PID profiles from these five different conditions in Figures 2A and S1. Model parameters are λ = 0.02, αtranscription = 2, αtranslation = 2, and ki = 50 (for all i = 1, 2, 3). Relative activation rates for transcription, αm, depended on the number of activating and inhibiting regulators present in each possible state Sm:αm = 0.1 (for the unbound state, i.e., basal transcription), 0.001 if an inhibitor was bound (we assumed inhibition dominated activation), and 5 if only activating regulators were bound.

60The second mass-action model that we consider also includes seven species: genes (gi) and mRNA/protein (xi) for i = 1, 2, 3 and the stimulating ligand s. We assume that protein and mRNA concentrations are equal (i.e., translation is instantaneous) and, unlike our first model, assume that a gene can only be bound by a single protein at any time. The possible reactions and associated propensities are

61for basal transcription and protein decay, protein binding and unbinding from a target gene, and transcription from a gene bound to a regulating protein, respectively (where gixj indicates gene i is bound by regulating protein j).

62For the mass-action model, we simulated time courses from time 0 to 400, with the stimulating ligand present from time 200 at a constant level of 20 molecules, and recorded the system state at 21 equally spaced time points. We repeated the simulations 50 times with the stimulating ligand targeting each gene to give a total of 1,050 data points (cells) that are used to calculate PID measures. We initiated simulations with two copies of each gene (gi = 2), no stimulating ligand, and initial mRNA/protein levels sampled from a uniform distribution (xi∼U(0, 50)); we perform simulations for five different initial conditions and calculate the mean PID profiles from these five conditions (plotted in Figures 2A and S1). Model parameters are ktxn = 1,kdecay = 0.05, kon = 0.01, koff = 0.25, and kregulated = 10 or 0.1 for activating and inhibiting regulation, respectively. Information measures were calculated from these data using the MATLAB package written by, following discretization with the AstroML implementation of the Bayesian blocks algorithm.

63Data are simulated using GeneNetWeaver, a software package that generates stochastic simulations from in silico networks that are designed to be representative of real biological network structures (they are created by extracting subnetworks from known E. coli and S. cerevisiae transcriptional networks). This software has become a common tool for simulating gene expression data — including its use as part of several DREAM (Dialogue on Reverse Engineering Assessment and Methods) network inference competitions — which aims to provide unbiased datasets that do not favor particular inference methods and networks that retain characteristics of real GRNs. We compare the network inference algorithms using ten networks, five with 50 genes and five with 100 genes; for each network size, there are two E. coli and three S. cerevisiae networks, with average node degrees ranging from 1.19 edges per node to 5.51 edges per node.

64GeneNetWeaver uses dynamical models that consider mRNA transcription and translation processes and generates time-series simulations using stochastic differential equations to model dynamical noise and a mixed normal and log-normal model to represent microarray noise. To mimic single-cell data, we simulated thousands of time-series experiments for each network, using the default settings, with mRNA measurements generated according to the default settings and the default time points: times 0 to 1,000, in steps of 50. We sampled a single time point from each time series, representing a single cell: for the large datasets, we sampled 100 cells from each of the time points (2,100 cells in total); for the medium datasets, we sampled 100 cells from time 0 onwards in time steps of 150 (700 cells); and for the small datasets, we sampled 20 cells from all time points between 0 and 300 (140 cells). Where the dataset size is unspecified, we have used the medium datasets. For our initial exploration of PID measures, we also used data simulated from a 10-gene network (shown in Figure 3A); in this case we again used the temporal sampling scheme for medium datasets.

65In order to test the robustness of our method to zero-inflated data, typical of single-cell experiments, we further simulated dropout datasets. Zero measurements appear to be a combination of technical errors and genuine lack of expression due to stochasticity or biological state and are more common in transcripts with low abundance. We simulated dropout events at two rates, such that expression values in the lowest 50% (high rate) or 20% (low rate) for each gene had a 50% probability of being recorded as 0.

66AUROC and AUPR curves are calculated by comparing the inferred networks (which assign a score to every potential network edge) with the true network used to simulate data and identifying the numbers of correctly (and incorrectly) assigned edges as the threshold for edge inclusion is varied. AUROC is calculated from the area under the ROC curve, which is a plot of the false-positive rate (FPR) on the x axis versus the true-positive rate (TPR) on the y axis. AUPR is the area under the curve for a plot of precision (y axis) versus recall (equal to TPR) on the x axis. These quantities are calculated aswhere TP and FP indicate the numbers of true and false positives, and TN and FN are true and false negatives. For networks where the number of negatives is much greater than the number of positives, AUPR is considered a better metric for comparing algorithm performance.

67We use these scores to compare the performance of our method relative to existing inference algorithms (Figures 5 and S3). We used the implementations in the R package minet for the existing inference algorithms; we used the default or suggested parameters within this package except for the tolerance, τ, for ARACNE for which we either used the default (for the results in Figure S3) or 0.1 (for the results in Figure 5) as recommended in the original publication describing this algorithm.

68We analyzed three published qRT-PCR datasets to illustrate our network inference algorithm. Normalized Ct values from were subtracted from the assumed maximum, 40, and the resulting dCt values (for 87 genes and 681 cells) were used in our analyses. Normalized dCt values from are used directly for our analyses; we used data from the 20 genes they represent in their network model and 3,934 cells. Raw Ct data from are treated as described by the original authors (dCt values are calculated assuming a limit of detection of 28 and normalized on a cell-wise basis by subtracting the mean expression of housekeeping genes Actb and Gapdh; all values corresponding to expression below the limit of detection are set to −15); we used data from 46 genes (i.e., we excluded the housekeeping genes used for normalization) and 442 cells.

69The InformationMeasures.jl package is available from https://github.com/Tchanders/InformationMeasures.jl.

70A Julia package for running the PIDC, PUC, CLR, and MI algorithms is available from https://github.com/Tchanders/NetworkInference.jl.

71Tutorials and simulated datasets are available from https://github.com/Tchanders/network_inference_tutorials.

72T.E.C., M.P.H.S., and A.C.B. designed and performed the research, analyzed data, and wrote the paper. All authors read and approved the final manuscript.