Abstract
We employ the language of Bayesian networks to systematically construct gene-regulation topologies from deep-sequencing single-nucleus RNA-Seq data for human neurons. From the perspective of the cell-state potential landscape, we identify attractors that correspond closely to different neuron subtypes. Attractors are also recovered for cell states from an independent data set confirming our models accurate description of global genetic regulations across differing cell types of the neocortex (not included in the training data). Our model recovers experimentally confirmed genetic regulations and community analysis reveals genetic associations in common pathways. Via a comprehensive scan of all theoretical three-gene perturbations of gene knockout and overexpression, we discover novel neuronal trans-differrentiation recipes (including perturbations of SATB2, GAD1, POU6F2 and ADARB2) for excitatory projection neuron and inhibitory interneuron subtypes.
The classification of cortical neurons is a debated topic with differing schemes using anatomical, molecular and physiological characteristics in order to make distinctions. It is generally accepted that there exist two major groups of neurons namely, excitatory Projection Neurons (PNs) (Greig et al. 2013) and inhibitory Interneurons (INs) (Kepecs and Fishell 2014). All neurons are generated only during embryonic development (Lodato et al. 2015) after which class-specific traits remain unchanged for the life of the organism. It is classically thought this precludes any change in identity postnatally. Intriguingly neurons may exhibit more plasticity than previously thought. As far back as 2002, astrocytes were directly reprogrammed into neurons and more recently post-mitotic neurons have been converted from one subtype to another in young animals as reviewed by Amamoto et al. (Amamoto and Arlotta 2014).
A Bayesian network (BN) is a graph-based model of joint multivariate probability distributions that captures properties of conditional independence between variables. Bayesian networks can be used for representing statistical dependencies in a set of data and were applied to the problem of reconstructing gene regulation networks (GRN) from expression data by Friedmann et al. (Friedman et al. 2000) and Hartemink et al. (Hartemink et al. 2001). It is known that the protein transcription factor produced by one gene can have a causal effect on the expression of another gene. BN can be used to represent the conditional dependencies between genes and thus interpret these as causal patterns of gene regulations.
We challenge the paradigm that neurons of the mammalian cortex are a permanently post-mitotic and differentiated cell type via modeling genetic perturbations that facilitate direct transdifferentiation. In this theoretical study we present the application of Bayesian network techniques to high quality deep-sequencing data in order to reverse engineer the genetic regulations in human neurons. We identify those attractors that correspond to different neuron subtypes and validate our model with an independent data set. Using dynamic bayesian inference we derive interconversion recipes between differing neuron subtypes and from the perspective of the cell-state potential energy landscape, describe those interconversion pathways.
Method
Data processing, clustering and discretization
Lake et al. (Lake et al. 2016) previously conducted single-nucleus RNA sequencing on post-mortem adult human cerebral cortex and generated 3,227 quality-filtered single neuron data sets. These nuclei were subsequently resolved into 17 clusters, based on the differential regulation of 16,242 protein-coding genes, through repeated rounds of unsupervised hierarchical clustering and supervised classification (technical details can be found in ref. (Lake et al. 2016)). Figure 1 A.i. is a reproduction of the authors hierarchical tree down to Level 2 including the clusters considered in this work. At each of the 3 splits, we consider the ten-fold differentially expressed genes (DEGs) giving a total of 74 unique genes. Transcription levels were previously analyzed as 
of transcript per million mapped reads (TPM). Thus, for the complete data set at Level 2 (
samples) we calculate the weighted arithmetic mean for each gene 
where 
where n is the number of samples in each cluster. Each data point was subsequently discretized according to:
(1)thus transforming the data:
where gene x runs from 1 to 
and sample j runs from 1 to 
The experimental barcodes for each cluster, post-discretization, are shown in Figure 2.
Workflow pipeline. A.ii. Hierarchical clustering of 3,227 quality-filtered single neuron data sets from previous single-nucleus RNA sequencing study. Number of 10-fold up-(↑) and down-regulated (↓) DEGs given at each junction. These are defined as up-regulated in the subsequent left hand branch and down-regulated in the left hand branch each relative to the right hand branch. Terminal clusters I - IV at Level 2 used in this work.(Lake et al. 2016). A.i. 77 ten-fold DEGs (of which 74 are unique) used to train networks. Splits refer to junctions in the hierachical tree in Figure 1 A.ii. Genes in bold occur in multiple splits. B. Data discretization on the weighted arithmetic mean of the log2(TPM) for each gene across 3,203 Level 2 single neurons. Data downsampled to 1,920 samples. C. Structure learning directed acyclic gene reulation networks using the discretized downsampled data with local and global optimization routines. 20 random seeds used to generate 20 different structures. 5 structures randomly chosen for further calculations. D.i. 450,000 random initialisations of the nodes in the continuous interval [0,1] for the 5 network structures. 2 TBN DBN inference performed for each network in each initial state. Converged attractor states subsequently clustered. D.ii. Nodes in 5 network structures initialised in the continuous interval [0,1] corresponding to the four neuronal subtype cell states. For each network structure in each initial state, 2 TBN DBN inference carried out for all 3-gene perturbation combinations (clamping nodes as overexpressed or knocked out for duration of inference). Subsequent node probabilities averaged over the 5 structures for each 3-gene perturbation in each initial state.
Experimental barcodes for clusters I - IV post data processing. For gene expression probabilities 
code displayed as green and for gene expression probabilities 
code displayed as red. Genes that display the same pattern across all four clusters are displayed on the same row. Genes are colored red, purple, green and blue for those expressed in clusters I - IV respectively. DEG that are expressed in both PN subtypes I and II are colored magenta and DEG expressed in both IN subtypes III and IV are colored cyan.
In order to prevent the network topologies being biased toward any one neuronal sub-type, the discretized data were down-sampled via the random removal of samples from clusters I, II and III until each cluster contained 480 samples, thus matching the lowest cluster size, that of cluster IV. For each of the three downsampled clusters, the appropriate number of samples were randomly selected for removal and the probability of expression (assigned according to the fraction of the 480 samples that were in a state 1 post-discretization) was calculated for each gene. Comparison was then made to the probability of expression for each gene from the complete, non-downsampled cluster via the root mean square deviation (RMSD) and Pearson R correlation summed over all 74 genes. This process was repeated 
times for each cluster and the set of samples that gave the lowest resultant RMSD were removed. The 
values for each downsampled cluster compared to the relevant complete cluster were 
and 
This gave a data matrix of size 
for 20 separate structure learning runs.
Structure learning
For a given directed acyclic graph (DAG) model G based on data D with n binomial variables, it can be shown that
(2)where n is the number of variables, 
is the number of instantiations of the parents 
is the ascertained prior belief of the number of times 
takes its first value when parents 
are in their 
instantiation, 
is equivalent to 
but with 
taking its second value, 
is the number of times in the data 
takes its first value when parents 
are in their 
instantiation, 
is equivalent to 
but with 
taking its second value, 
and 
is the gamma function (Neapolitan 2009). Equation 2 is defined as the Bayesian score assuming Dirichlet priors. In order to punish overly complex DAGs and reduce the possibility of overfitting, we use the Bayesian information criterion (
) to score structures:
(3)which includes an error term, where m is the number of samples and d is the dimension of the DAG i.e., the number of parameters.
Our procedure follows a three-stage score-based approach common to the method of Chang et al. (Chang et al. 2011) but omitting prior knowledge incorporation and is thus purely data driven. We give a brief overview here for the readers convenience. Since the problem of learning optimal structure is NP-hard (Chickering et al. 1995) we use heuristics in the form of the greedy algorithm to maximimise the 
score during the first stage. Starting with an empty network, two random nodes (genes) A and B are selected. If no edge exists between them, either the directed edge A → B or the opposite regulation B → A is generated, each with a probability of 0.5. If an edge already exists between the two nodes (a possibility from the second step onwards), it is either reversed or deleted each with a 0.5 probability. For each of the four outcomes, the change to the network is accepted if the 
score increases, else it is rejected. In our case, this started with an empty network and was iterated 
times (
NSteps).
During the second stage, we employ the metaheuristic approach of simulated annealing in order to approximate the global minimum. The network is instantaneously heated to a temperature T and uniformly cooled over the course of the stage. According to the logic set out in stage 1, an edge is generated. However, if this leads to an unfavorable decrease in the 
score, the edge is accepted with a probability of 
Figure 3 shows this probability as a function of the 
for a given edge introduction, at different temperatures. Stage 2 was iterated 
times (NSteps) from a starting 
down to 
Acceptance criterion for edges that are unfavorable to 
score during simulated annealing stage of structure learning. Representative temperatures in the range 
plotted.
Stage 3 consists of a final stage 
maximization via the greedy algorithm, identical to the protocol of stage 1. This was iterated 
times (
NSteps). The change in the average 
score was converged with the value of NSteps as shown in Figure 4.
Convergence of 
(referenced from structure with 1.1M Nsteps) with NSteps (defined as number of steps during simulated annealing phase). Each value is the arithmetic mean over 5 independent structure runs.
Since the data used was non-temporal, only edges that did not introduce a loop into the structure were accepted. We employed this structure learning procedure using 20 random seeds and learnt an ensemble of models.
Inference
During inference we apply Bayes’ rules to obtain the posterior probabilities. Since we are interested in developing a gene regulation model that can be initialised into neuronal subtype cell states; for each source cluster I - IV (see Figure 1 A.i.), initial probabilities for each node x were defined as the probability that the node was in a state 1, i.e., expressed, for the given cluster and were input as values in the continuous interval 
These probabilities were assigned according to the fraction of the total samples, in the given cluster, that where in a state 1 post-discretization.
Model averaging
We use the full Bayesian framework and do not attempt to approximate one true underlying distribution with a single structure. The a posteriori distribution of models is:
(4)where 
is the likelihood of the model and the prior, 
is not assumed to be uniformly distributed and thus not constant. We therefore perform probabilistic inference by model averaging where the averaged conditional distribution of variable X is obtained by integrating over models:
(5)Since the model space is written as 
where 
is the continuous ensemble of the conditional probability table (CPT) configurations for each structure 
For every 
each possible parameterization in the CPT configuration ensemble 
defines a member, 
and the distribution is normalized against all models:
(6)where the normalization factor 
Equation 5 is thus extended to a double-integral over structure space and structure-dependent parameter space:
(7)The Bayesian network models k sum over discrete structure space and the parameter vector configuration θ integrates over continuous parameter space and can thus become intractable by analytical methods. Within the Markov Chain Monte Carlo (MCMC) approach we use the order statistics of a uniform distribution [0,1] to simulate a sample from posterior distribution 
(a beta-density function).
Learned parameter retrieval
For a given structure, the probability of gene X with no parents remains as the initial probability for the course of the inference run. For gene X with one parent gene Y, the probability of expression of gene X (i.e. 
) is dependent on the two parameters 
and 
which represent the two respective conditional probabilities. We calculate these conditional probabilities according to the associated beta-density function as:
(8)and
(9)where 
are the number of data samples where Y and X occur with the values h and k in the discretization 
The additive value of 1 is for the case that there are no data samples in the given instantiation. A gene X with multiple parents 
has 
associated with it, where 
is the number of parents of gene X. Thus the more general formulation for the parameters is:
(10)where 
is the number of samples in the data where the parents of X occur collectively in state j and X occurs in state k 
Dynamic Bayesian network
To simulate the evolution of cell state from initial state to the equilibrium or steady state solution we use the dynamic Bayesian network (DBN) model. Here the probabilistic inference is performed by using a 2-time slice Bayesian network (2TBN) and the interface algorithm, (Murphy 2002) which uses static junction trees as a subroutine to compute exact inference in the 2TBN which is then repeated sequentially over time. Accordingly node (gene) probabilities (expressions) evolve over time according to:
(11)where 
is the product of the probabilities of each parent of X being in the binary state 
corresponding to the collective state j. Thus the probability of X being expressed at time t depends on its parameters, 
as well as the probabilities of expression of its parents 
at time 
Each set of parameter samples θ forms an instance of the DBN model which were averaged for each inference run. The inference runs for each structure 
were then averaged for the final result. During perturbations, node combinations were either clamped at 
(overexpression) or 
(knockout) for the duration of the inference run.
Landscape analysis
In order to identify those attractors in the unperturbed state, we randomly initialise the genes in the continuous interval 
and perform 2TBN DBN. The resulting method corresponds to converging on the minima accessible to the region of cell state space in which the initialisation took place. As per the method of inference under perturbation, node probabilities evolve over time according to Equation 11. In the case of attractor identification we do not apply any clamps and all nodes are free to evolve to a steady state solution. For all inference calculations a set of parameter samples θ forms an instance which are averaged over for each unique random initialisation. This process was repeated 450,000 times until the number of unique attractors had converged (234 attractors defined in binary cell state space) thus the sample is large enough to accurately describe the entire space. The basin size of the potential energy landscape corresponding to each attractor can then be calculated as a percentage of the total number of samples that converge to each state.
Transition states
Transition states during inference were calculated as per Chang et al. (Chang et al. 2011) by applying a maximum-a-posterior (MAP) estimation to predict the state-transition pathways. That is, at each time step t, the state which maximizes the cell state posterior at the current time step is selected as the current cell state:
(12)where the probability propagation in DBN cell states is defined as:
(13)For each unique binarised cell state, the state probability for the ith state (
) and thus the potential energy, 
of the state i can be calculated. These potential energy values for all 
binary states can be represented as 
Cell potentials differ according to cell state conditions. In this work we investigate the cell state potential changes for specific transition states under a specific 3-gene perturbation 
Method validation
Buganim et al. (Buganim et al. 2012) have previously conducted a single-cell gene-expression analysis of mouse embryonic fibroblasts (MEFs) during cellular reprogramming. They profiled 48 genes from early time points, intermediate cells, and fully re-programmed iPSCs. These data were used to train a simplified Bayes model of hierarchical gene regulation in iPSCs. Using their regulation topology they chose five transcription factor combinations predicted to induce activation of the pluripotency circuitry and generate fully reprogrammed iPSCs. These were experimentally verified via flow cytomeric analysis using OCT4-GFP with 
reprogramming efficiency. Using this independent dataset we learnt 20 independent structures and performed inference calculations to predict the effect of their experimentally verified three and four gene overexpression combinations using our methods. The reference Pearson R correlation between the MEF initial cell-state and the fully reprogrammed iPSc final state was calculated to be 0.155. The average Pearson R correlation of the experimentally verified reprogramming recipes and the final iPSc cell-state was predicted to be 0.838 using our methods. This was compared to 20 random 3-gene overexpression combinations with an average Pearson R correlation of 0.270.
Data availability
Data used has been previously deposited with dbGaP (accession phs000833.v3.p1), curated by the NIH Single Cell Analysis Program Transcriptome (SCAP-T) Project (http://www.scap-t.org) as stated in Lake et al.(Lake et al. 2016). Supplemental material available at Figshare: https://doi.org/10.25387/g3.6349553.
Results And Discussion
Experimental expression profiles
Lake et al. (Lake et al. 2016) categorized excitatory PN by layer position and as such, cluster I neurons (which were further split in their hierachical tree) were labeled as a combination of granular neurons (GN) from layer 4, sub-cortical projection neurons (SCPNs) from layer 5 and cortico-thalamic projection neurons (CThPNs) from layer 6. Cluster II were classified as cortical projection neurons (CPNs) residing in layers 2/3. The CPNs were shown to express CPN-associated CUX2 and the layer 2/3 marker gene LAMP5 both of which were 10-fold DEGs between clusters I and II and thus included as nodes in our network. Functionally, layers 1-3, termed the supragranular layers, are unique in the neocortex and are the primary origin and termination of intracortical connections. These can be functionally contrasted with internal granular layer 4 and infragranular layers 5 and 6. RORB a marker for layer 4 neurons (Schaeren-Wiemers et al. 1997) is also shown to be up-regulated in cluster I compared with cluster II consistent with this analysis. Interneuron subcategories were found to be distributed across across the neocortex and were classified based on developmental origin. (Lake et al. 2016) Cluster IV IN were found to originate from lateral (LGE), or caudal ganglionic eminences (CGE) and were VIP+ and RELN+ with positive expression of P8 and NR2F2. Whereas cluster III IN showed MGE marker expression such as LHX6 and SATB1.
The 74 unique genes used to train the GRNs in this work are given in Figure 1 A.ii. They are those that are 10-fold DEGs between the clusters at Levels 1 and 2 of the hierarchical tree (see Figure 1 A.i.). The retrieval of expression profiles for clusters I - IV post data processing and discretization was conducted via the generation of experimental barcodes (see Figure 2). For each cluster, the gene probabilities were calculated as outlined in the Methods section and subsequently binarised based on a cutoff of 0.5. In this way each cluster is represented as one of the 
possible states. As a further assessment of the clustering and cluster uniqueness, we compared the RMSD and Pearson R correlation, summed over all node probabilities in 
between all four clusters (see Table I). From this we retrieve the fact that inter-neuron expression differences for the excitatory PN clusters I and II are less distinct (
and 
) than those between the IN clusters III and IV (
and 
) in agreement with Lake et al. (Lake et al. 2016). For all other intra-neuron subtype comparisons (i.e., PN sub-categories vs. IN sub-categories) we find negative correlations exist in the range 
to 
thus find, with the possible exception of comparison between clusters I and II, the cluster expression profiles to be adequately unique (within our state space defined by the 74 DEGs) after data processing and discretization.
Landscape analysis
Due to the computational intensity of 2TBN DBN inference over a sufficiently large sample of state space, five structures were randomly chosen, from the ensemble of 20 BN structures that were learnt, for landscape analysis and inference calculations. For attractor analysis node probabilities for the 74 genes (1 single cell state) were randomly initialised in the continuous interval 
and the node probabilities were converged to a steady state. Post ensemble averaging the node probabilities were discretized. The process was repeated for more randomly selected initial cell states until the number of unique attractors (local minima) was converged. We found that the random sampling of 450,000 initial states was sufficient to identify all the major attractors in the network. 1082 unique attractors were found and hierarchically clustered using the heatmap.2 function with the Euclidean distance metric in R as shown in Figure 5. Arbitrarily cutting the tree at a distance of 1.3, groups the attractors into 3 representative cell state clusters. The basin size of a given attractor on the cell state potential landscape can be defined as the percentage of random initial cell states that converge to the given attractor. The basin size for each unique attractor within each of the three representative cell state clusters A, B or C was summed to give the total basin size for those representative state clusters.
Attractor heat map for all 1082 unique attractors of the unperturbed network model. Cell states defined using quaternary intervals 
and hierarchically clustered using the heatmap.2 function in R. Basin sizes for the 3 representative cell state clusters A, B and C given in legend.
An important test of the four models was to assess the correlation between the initial experimental gene expression values (calculated as described in the Method Section) and those expression values (node probabilities) to which the structure relaxes during 2TBN DBN inference in the unperturbed state. In terms of the cell-state potential landscape, this represents the proximity of the nearest (defined in terms of cell state similarity) “local” minimum or attractor to which the inference results converge. All Pearson R rank correlation coefficients between initial experimental values and relaxed unperturbed inference results are given in Table II. All values are in the range 
showing that the model adequately describes the neuronal cell states corresponding to subcategories I–IV.
The two unique attractors corresponding to the two PN subcategories I and II cluster together in representative cell state cluster C (see Figure 5) with a basin size covering 27% of cell state space. Further, the two unique attractors corresponding to the two IN subcategories III and IV cluster together in representative cell state cluster B with a basin size covering 28% of cell state space.
In order to further validate our network model, additional single-cell RNA-Seq samples, taken from further Brodmann areas in the adult human cerebral cortex and processed in identical fashion to Lake et al. (Lake et al. 2016), were used. These samples included neuronal and non-neuronal single cells, such as glia, astrocytes, oligodendrocytes and microglia. Subsequent sample filtering resulted in 546 quality samples which were discretized in the binary interval [0,1] based on the arithmetic mean for each gene summed over all samples. This resulted in all 546 samples displaying a unique cell state (using the same 74 gene space as the model). To make comparison to the attractors predicted by the network model, the 1082 quaternary cell states (clustered in Figure 5) were binarised resulting in 234 unique attractors. Comparison between the additional samples and network attractors was made with the Hamming distance metric. It was found that 451 of the 546 additional samples (83%) had a hamming distance of 20 or less with one or more of the 234 network attractors, i.e., 54 or more genes (73%+) were in the same state of expression. This is suggestive of the fact that the network model (trained only on neuronal subtypes) captures global genetic regulations for other non-neuronal cell types of the mammalian cerebral cortex in addition to accurately describing neuronal cell states.
Topological analysis
As previously described, 20 BN structures were trained from the combined down-sampled data of clusters I - IV (see Figure 1). The mean number of edges learnt across all 20 structures was 
and the standard deviation 
The relative sparsity of these networks (
) owes to the inclusion of the error term in the 
scoring function (see Equation 3) leading to the penalisation of overly complex structures. (We use, in the standard definition of network density (D), the number of possible edges in a complete graph (Emax = n(n − 1)) in a directed network where we do not allow self-regulation/loops but hypothetically cycles would be included (a formulation forbidden in our BN learning approach).) A consequence of this and learning multiple structures is that despite the highly stochastic nature of the learning protocol, edges that do occur in higher frequencies across structure learning runs should represent true biological regulations that are coded for in the data.
Community Detection & Pathway Analysis
Figure 6 shows the merged GRN for all 20 structures with only those edges that occur in 
of structures displayed. Node sizes are scaled by out-degree. Community detection analysis using the fast unfolding heuristic algorithm of Blondel et al. (Blondel et al. 2008) shows there to be 4 communities at a resolution of 1.55. Interestingly of the 20 nodes depicted as being members of the orange community 13 (65%) were also up-regulated IN relative to PN and as such were DEGs from the root split, furthermore 6 (30%) of genes in this community were further up-regulated in cluster III relative to cluster IV from the IN split. This suggests that regulatory mechanisms in this community likely lead to strong co-expression and activating regulation between the nodes in IN and more specifically cluster III IN (MGE derived). Furthermore 13 of the 16 genes (81%) in the purple community are up-regulated in cluster IV IN (LGE/CGE derived) relative to cluster III IN. The same analysis on the green community reveals that 12/23 (52%) genes are broadly up-regulated in PN at the root split and 6/23 (26%) are further up-regulated in cluster II CPN. Finally, the blue community is almost exclusively, with the exception of two nodes, from the PN split with 5 and 7 genes up- and down-regulated respectively in the I GN/SCPN/CThPN cluster. Based on the fact that the network gene list derives from IN and PN specific DEGs, the neuronal subtype DEG specificity, as related to the communities, is partially expected.
Merged GRN for 20 independent BN structures. 228 unique edges with frequency 
included (from a total of 870 unique edges summed over all 20 BN). Displayed using Yifan Hu algorithm as implemented in Gephi version 0.9.1. Node size proportional to out-degree and edge width proportional to frequency. Community detection algorithm (Blondel et al. 2008) run with resolution of 1.55. Modularity = 0.421. Number of communities = 4.
Pathway analysis of all 74 genes reveals a significant number of enriched pathways from “Signal Transduction””, “Immune System”, “Transmembrane transport of small molecules” and “Neuronal System” among others. In particular within the “Neuronal System”, the “Neurotransmitter release cycle” (P value 
) includes the genes GAD1, GAD2, SLC6A1 and SLC17A7, the first 3 of which are part of the orange community and form a clique (complete subnetwork with edges between all 3 members in both directions) each forming an edge with each other at a frequency of 
across all 20 structures. The pathway “Signaling by Type 1 Insulin-like Growth Factor 1 Receptor” (P value 
) containing the genes IGF1, ERBB4, KIT and EGFR is also enriched. All of these genes are identified by community analysis to be part of the same community and further, the two edges 
and 
occur in 
and 
of all structures respectively.
Edge distribution
The top 27 edges that occur in all 20 structures learnt are given in Table III. One of these regulations is 
SLC17A7, a known regulator of brain physiology, is a brain-specific solute carrier and is found to regulate the neuronal signal transducer, chimerin-1 (CHN1). SLC17A7 specifically functions as as glutamate transporter and it has been found that 
-chimerin regulates dendritic spine density. (Van de Ven 2005) Spine morphological changes, associated with long-term depression, can be induced in hippocampal neurons by metabotropic glutamate receptor activity suggesting possible support for this learnt regulation.
The POU family members are transcriptional regulators, many of which are known to control cell type-specific differentiation pathways. STXBP6 codes for the syntaxin binding protein 6 (amisyn) and an edge that occurs in all structures is 
STXBP6 is known to regulate SNARE complex assembly, a protein complex involved in membrane fusion, that play an important role in neurotransmitter release.
ADARB2 forms two edges that occur in all structures and is the rank 3 node as ranked by total degree (see Table IV) suggesting it plays an important role in neuronal gene regulation and identity. It is a member of the ADAR family which contains 3 members, two of which, ADAR and ADAR1, are catalytically active. The ADARB2 gene encodes a catalytically inactive protein, expressed in brain, amygdala and thalamus. (Hogg et al. 2011) It is known to prevent the binding of other ADAR enzymes to targets in vitro, and decreases the efficiency of these enzymes. These enzymes are responsible for RNA editing via the conversion of adenosine to idenosine which has been observed in some pri-miRNAs(Kawahara et al. 2007; Yang et al. 2006); that can in turn affect the function of miRNAs which are thought to have a functional roles in gene regulation. The edge between 
occurs in all structures. Interestingly EGFR/MAPK has been shown to regulate AGO2, (Adams et al. 2009) which itself is a member of the AGO protein family that play a central role in the function of the RNA-induced silencing complex (RISC) and therefore potentially miRNA function. (Höck et al. 2007)
NFIB, a gene essential for brain development in mice (Steele-Perkins et al. 2005) and NFIX form another high frequency edge 
(see Table III). Both genes belong to the NFI family encoding site-specific transcription factors whose functional diversity is generated in part through protein heterodimerization, (Liu et al. 1997) thus providing strong evidence for a protein-protein interaction and a mechanism of co-regulation.
The edge 
occurs in 
of structures and is part of the purple community found to be up-regulated in cluster IV LGE/CGE derived IN. Putative homologs of these genes were found interacting in other organisms such as the protein-protein binding interaction in Drosophila melanogaster of CG31692 and ics and in saccharomyces cerevisiae the protein-protein interaction between RAS2 and CYR1.
Transdifferentiation gene recipes
Throughout this section we refer to the “source state” as node probabilities that are initialised to the given probability of expression for the source cluster I–IV and the “target state” as the node probabilities of the final or target cluster I–IV.
Perturbations were applied as either overexpression (clamping the node probability to 1 for the course of inference) or knockout (clamping the node probability to 0 for the course of inference). A full scan of the three-gene recipe combinatorial space was conducted. The calculations were performed using five randomly selected structures from the 20 trained and final node probabilities were averaged over these structures under each perturbation. Three-gene recipes for the 12 interconversions were ranked based on the RMSD between all node probabilities for non-perturbed genes (71) in the perturbed source state post relaxation and the corresponding node probabilities of the target state. The five best recipes for each of the 12 interconversions are given in Supplementary Tables I and II.
Table V shows the best 12 interconversion recipes. We find symmetries exist between the recipes, for example with the exception of S-II→T-I, for conversion to PN subtypes I and II the overexpression (↑) of SATB2 is in all recipes (irrespective of source cluster type). Contrastingly, all of the best 6 conversion recipes to IN subtypes III and IV include the knockout (↓) of SATB2. SATB2 is a DNA-binding protein that regulates chromatin organization and gene expression and is important in the development of corticocortical connectivity in the developing cerebral cortex in mice. (Alcamo et al. 2008) Broadly defined as an excitatory marker, SATB2 was found to be regulated between the PN and IN and is a DEG at the root split (see Figure 1 A.ii.). Topologically, SATB2 is the rank-two gene by degree (see Table IV) forming 18.15 edges on average per structure and further occupies a central position in the network making connections to all 4 communities (see Figure 6). The SATB2 protein through its interactions with both the CTIP2 promoter upstream region and histone deacetylase complex, controls chromatin remodeling. Upper layer pyramidal neurons lose their identity in the absence of SATB2 (Britanova et al. 2008) perhaps consistent with our prediction of SATB2 knockout being important in inter-neuron PN → IN transdifferentiation and also a regulation that is required to be “held in place” for inter-neuron IN subcategory → IN subcategory transdifferentiation. Similarly, but in reverse regulation logic to SATB2, the knockout of GAD1 is in all but one of the 6 recipes for conversion to PN subtypes and its overexpression is in 50% of the recipes for conversion to IN subtypes namely, S-I→T-III, S-II→T-III and S-II→T-IV. GAD1 is an inhibitory marker and is up-regulated in IN. It is the rank one node by degree forming 18.20 edges on average per structure and connects to nodes from three of the four communities identified (see Figure 6). GAD1 is involved in pathways including “Neurotransmitter release cycle” and “Transmission across Chemical Synapses” and is an integral enzyme in “Gaba Synthesis”. The overexpression of GAD1 is consistent with transdifferentiation to IN targets since the majority of IN are GABAergic.
In addition to these inter-neuron transdifferentiation recipe symmetries there also exist symmetries in targeting specific PN or IN subtypes. For example, the best 3 recipes targeting T-I include the overexpression of POU6F2 and by contrast all the best 3 recipes to T-II include the knockout of POU6F2. POU6F2 is a transcription factor involved in DNA binding and is only expressed in the CNS. Moreover the gene is enriched in GO terms for “central nervous system development” and “regulation of transcription, DNA-templated” consistent with its high rank by degree and the fact that it forms edges as a parent node to nodes in three communities. Finally of note in the best transdifferentiation recipes (see Table V) is ADARB2, of which the knockout and overexpression targets cluster III and cluster IV IN respectively. This is the rank three node by degree and is the most connected gene in the purple community in Figure 6 and its functions were discussed in the Edge Distribution subsection.
Figure 7 shows representative plots for the three best 3-gene recipes from source cluster S-I. Each scatter plot contains the target state node probabilities vs. first the source state node probabilities in the unperturbed (experimental vs. experimental) states in cyan and second the perturbed (target experimental vs. source theoretical perturbed) states in dark blue. We can see remarkable reprogramming success as reflected by the improvement in Pearson R correlation coefficient which in the case of S-I→T-IV changes from a strong negative correlation in the experimental vs. experimental plot of 
to 0.9589 under the perturbation ERBB4↑ / ADARB2↑ / SATB2↓. All the 12 best transdifferentiation recipes achieve final correlations in the range 
(see Supplementary Tables I and II).
The 3 best transdifferentiaion recipes as ranked by RMSD between source cluster I perturbed states and experimental target states. Subplots are target state node probabilities vs. source state node probabilities in the unperturbed (experimental vs. experimental) and perturbed (target experimental vs. source theoretical perturbed) states, represented by cyan and dark blue data points respectively. Best recipes are highlighted in green (overexpression clamp) and red (knockout clamp) and arrows are drawn to show the direction and magnitude of node probability for the given clamp.
Transition state analysis
It is instructive to monitor cell state probabilities during the transdifferentiation procedure via node clamping. In terms of the cell state potential energy landscape, each of the 
cell states (as represented by a unique combination of 74 binarised node states) in the unperturbed landscape has a potential energy associated with it that is calculated as 
where (
) is the state probability for the ith state. Under a 3-gene perturbation the available cell states with a finite probability is reduced by a factor 1/8 to 
and the probability of remaining states also changes. This adjustment to the landscape results in the re-positioning of minima and of energy barriers on the landscape which effectively makes states which were previously inaccessible, open to sampling. Figure 8 shows the probability of the states along the transition paths for 10 independent runs for the transdifferentiation S-I→T-IV under the perturbation ERBB4↑ / ADARB2↑ / SATB2↓ for one of the five structures used in inference. The initial unperturbed state (labeled transition state 1) has a potential energy of 33, the perturbation is then applied which raises the energy of the system to 
“transporting” the cell to a new, previously inaccessible area on the potential energy landscape. This new energy allows the cell to now relax to the new minima which coincides with the target cluster T-IV. We can see that the probability is converged in 
transition states which is common for most recipes.
Cell state potential changes along a representative path for a single structure under perturbed conditions. Each line represents a different sampling of the conditional probabilities. The best 3-gene recipe for the interconversion between S-I→ T-IV, ERBB4↑ / ADARB2↑ / SATB2↓, is shown.
Conclusions
In this work, we applied Bayesian network methods to make de novo predictions for neuronal transdifferentiation recipes between Projection neuron and Interneuron subtypes. Our network, trained on high quality single-cell RNA-Seq data, accurately describes the four cell subtypes in the unperturbed state and is well validated against an additional data set of single-cells from more varied areas of the human cerebral cortex, using attractor analysis. Many of the regulatory edges learnt between the genes are validated from the wider literature and community analysis reveals significant enrichment in neuron specific pathways among others. We conducted a systematic search for transdifferentiation recipes that could achieve reprogramming. The three-gene recipes identified achieved remarkable success the best of which achieve final correlations, with the target state, in the range 
Master inter-neuron regulators are identified as SATB2 and GAD1 as well the identification of POU6F2 and ADARB2 as important intra-neuron regulators.
Acknowledgments
The authors would like to thank Dr Richard Stein for the development of the structure learning and DBN inference code which was modified in this work.
Footnotes
Supplemental material available at Figshare: https://doi.org/10.25387/g3.6349553.
Communicating editor: C. Myers
- Received March 12, 2018.
- Accepted May 21, 2018.
- Copyright © 2018 Ainsworth et al.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
