The critical role of composition-independent temporal and autism pathologic changes in the human brain transcriptome

The functions of human brains highly depend on the precise temporal regulation of gene expression, and substantial transcriptome changes across lifespan have been observed. The substantial transcriptome alteration in neural disorders like autism has also been observed and is thought to be important for the pathology. While the cell type composition is known to be variable in brains, it remains unclear how it contributes to the temporal and pathological transcriptome changes in brains. Here, we applied the transcriptome deconvolution procedure to the age series RNA-seq data of healthy and autism samples, to quantify the contribution of cell type composition in shaping the temporal and autism pathological transcriptome in human brains. We estimated that composition change was the primary factor of both transcriptome changes. On the other hand, genes with substantial composition-independent expression changes were also observed in both cases. Those temporal and autism pathological composition-independent changes, many of which are related to synaptic functions, indicate the important intracellular regulatory changes in human brains in both processes.


Introduction
The development and aging of human brains are complex processes, which are shaped by anatomical and molecular changes 1, 2, 3, 4 . Autism spectrum disorder (ASD), as a common neurodevelopmental disorder, disrupts the critical developmental processes and results in the disruption of cognitive functions 5,6,7 . With the emergence of high-throughput measurement of different molecules, dozens of studies have been conducted to characterize age-related molecular changes in human brains, especially at the transcriptome level 8,9,10,11 . Meanwhile, several studies have also been conducted in order to investigate the transcriptome alteration in ASD 12,13,14 .
The human brain, however, is a highly complex and heterogeneous organ comprised of numerous different cell types, including neurons and multiple classes of non-neuronal glial cells -such as astrocytes, oligodendrocytes, oligodendrocyte precursor cells and microglia -as well as vascular, such as endothelial, cells. Each of those cell types expresses a distinct set of genes 15 and plays a unique and essential role in the development and functions of the brain 16 . Different cell types are also known to show different spatial-temporal distributions. Such complexity thus raises the unanswered questions: how much of the age-related molecular change in human brains, or specifically the age-related transcriptome change, is the direct consequence of the cell type composition change? If composition changes contribute to age-related expression changes, what's the biological meaning of the rest? Similar questions about the molecular changes in neural disorders, including ASD, also remain opened.
In order to comprehensively answer these questions, an accurate estimation of cell type composition in healthy and autistic human brains is required. Although experiments including stained cell counting 17,18,19,20 and large scale single-cell RNA-seq 21,22 have the potential to provide these data, these methods are either too labor-intensive or costly at present. Thus, the computational method of inversing sample heterogeneity, i.e. deconvolution, is one of the best alternative solutions to estimate the mixing percentage of different cell types 23,24,25 . The application of transcriptome deconvolution to human brain was long limited by the lack of comprehensive transcriptome data for the varied cell types in human brains. The recent emergence of the human brain single cell RNA-seq data, covering all the main cell types in the human brain 21 , now renders this approach more feasible. and root-mean-square deviation (RMSD), both of which were similarity measurements, CIBERSORT out-performed among the three algorithms, while DR performed slightly better than QP (Supplementary Figure S2). On the other hand, the discrimination measurement (see Methods), suggested that the cell type gene expression identities were reconstructed in the most discriminating way by DR, while CIBERSORT performed better than QP (Supplementary Figure S2).
In addition to the postnatal brain transcriptome, we further applied the three algorithms to the human embryonic developmental brain RNA-seq data obtained from BrainSpan database (Figure 1 and Supplementary Figure S1). DR-estimated composition pattern successfully recapitulated the domination of fetal replicating neurons dominate before 12 post conception weeks (pcw) and its dramatic elimination after that. This observation, coupling with the increase of fetal quiescent neurons, well matched with the neuronal proliferation that occurred during four pcw to 12 pcw 30 . Moreover, the estimated compositions by DR algorithm, especially those of fetal and adult neurons, presented a successive pattern with the postnatal composition changes estimated above, which further indicated the reliability and robustness of the composition estimation. On the other hand, CIBERSORT, although out-performed both DR and QP deconvolution according to the similarity measurement, failed to reproduce those scenarios.

Composition change is the primary power shaping the human brain postnatal temporal transcriptome
To comprehensively study the contribution of cell type composition changes in the human brain temporal transcriptome, we dissected the bulk gene expression into the composition-dependent and composition-independent components, based on the composition pattern estimated by DR algorithm (see Methods). Their contributions to the overall gene expression variance were estimated. On average, the composition dependent component, which represents the upper limit of transcriptome variance caused by variation of cell type composition in human brains, explained 26.4% of the total variance in the human postnatal age series data ( Figure 2a). Focusing on the 5,119 genes with age-related expression (referred as age-related expressed genes, age test BH-corrected FDR<0.05), the composition-dependent component contributed 34.7% of the total variance which was significantly higher than the other genes ( Figure 2b, permutation test, P<0.001).
To further investigate the roles of the composition-dependent and independent component in shaping the temporal gene expression pattern, we calculated the variances explained by age (age-explained variance) for each of the two components separately. Interestingly, the relative contribution of age-explained variance from the composition component (55.8%) was much larger than the proportion of composition variance among total variance (34.7%) for age-related expressed genes, while the difference was much smaller for the non-age-related expressed genes (27.8% vs. 22.1%) (Figure 2b). Altogether, these observations implied that the change of the composition-dependent component, i.e. the cell-type composition changes, was the main power shaping the observed temporal expression in human brains.
Additionally, the proportion of age-explained variance in the composition-dependent component was dramatically higher for the age-related expressed genes than for other genes (median=66.4%, permutation test, P<0.001, Figure 2c); meanwhile for the same genes, a moderate but significant increase of age-related variance proportion was observed for the composition-independent component (median=28.9%, permutation test, P<0.001, Figure 2c). These observations not only point to the significance of composition-dependent changes, but also to that of the composition-independent changes, some of which may represent the changes of molecular features in one or several cell types, also participate in shaping the temporal transcriptome in human brains. It is worth to mention, that all the above observations were reproducible based on CIBERSORT-based composition pattern (Supplementary Figure S3), suggesting that it is not artifact due to the DR-based deconvolution procedure.

The age-related changes in the composition-independent components are important and well regulated
To better understand the biological significance of the age-related changes in the composition-dependent and composition-independent component explicitly, we applied age tests to each of the two components (see Methods), to identify genes with significant age-related changes in either component. 8,156 and 1,455 genes were found with age-related changes in the composition-dependent and independent component, respectively (Figure 3a, Supplementary Table S2). Both of the two gene sets were largely overlapped with the age-related expressed genes (Fisher's exact test, P<0.0001). However, no significant overlap was observed between them (Fisher's exact test, P=0.216, odds ratio=0.932), implying the independent contribution of the two components to the temporal transcriptome in human brains. CIBERSORT-based composition pattern provided similar estimation, resulting in 6,671 and 1,129 genes with age-related changes in the composition-dependent and independent component, respectively. Majority of those genes were overlapping with genes identified based on DR-based estimation (Supplementary Figure S4).
We further grouped the 8,719 genes with age-related changes in at least one of the two components into three categories: G1 -genes with age-related changes in both components; G2 -genes with age-related changes only in the composition-dependent components; and G3 -genes with age-related changes only in the composition-independent components. Interestingly, these three groups of genes showed distinct temporal expression patterns (Figure 3b). G1 and G2 genes, and especially the latter gene, showed higher expression levels in early postnatal development stages, while G3 genes were highly expressed in adult stages. Different groups of genes were also enriched in different cell types: G1 -adult neuron, G2fetal quiescent neuron, fetal replicating neuron and adult neuron, and G3 -astrocyte and endothelial cells (Figure 3b).
More importantly, the three groups of genes showed distinctive functional enrichments (DAVID 31 , Supplementary Table S3). In brief, G1 genes were enriched for synapse and translation-related functions. G2 genes were enriched in transcription regulation, protein degradation, and cell cycle. Lastly, G3 genes were significantly involved in extracellular regions and metabolism. These results indicated the distinct biological significance of the age-related composition dependent and independent changes.
The expression of these three groups of genes should have been modulated by certain regulatory mechanisms, such as transcription factors (TFs). To test this, we estimated the enrichment of TF binding motifs in the promoter regions of genes in each category. We observed significant excess of enriched TF binding motifs (hypergeometric test, P<0.1) in the G2 and G3 genes (permutation test, P<0.001) (Figure 3c). In addition, the expression of representative TFs of the enriched TF binding motifs in groups showed significantly better correlation (Wilcoxon test, P<0.05) with their targets in the respective groups than expected by chance (permutation test, P<0.001) (Figure 3c). TFs with TF binding motifs enriched in G2 genes, e.g. CUX1 and E2F1, were mostly negatively correlated with their targets and had been shown to be relevant to cell migration, cell cycles and neuronal development and maturation 32, 33, 34 . On the other hand, most of the TFs with binding motifs enriched in G3 genes, e.g. SMAD3, SREBF1 and NR2F2, were positively correlated with their G3 target genes, many of which had been reported participating in signal transduction and metabolism of astrocytes 35,36,37 .

Significance of composition-independent component in ASD pathogenesis
Autism spectrum disorder (ASD) is a common neurodevelopmental disorder. Previous transcriptomic studies have suggested dramatic changes of gene expression 12, 13, 14 . For instance, using the same human PFC postnatal age series RNA-seq data set as described above, coupled with RNA-seq data of 34 PFC samples of autistic patients, Liu et al. identified 1,775 genes with differential expression (DE) between ASD cases and unaffected controls 14 . It is worthwhile to estimate the contribution of cell type composition as well as to identify the composition-independent alteration which may represent the active regulatory alternation in ASD patients.
We applied CIBERSORT, QP and DR algorithms to estimate composition patterns in the autistic samples. The three algorithms provided similar results, and interestingly, all showed significant smaller proportion of neurons as well as larger proportions of non-neuron cells (Figure 4a and Supplementary Figure S5). The amplitudes of composition changes were significantly correlated among the estimations based on the three algorithms (Supplementary Figure S5), eliminating the possibility of algorithm artifact. The increase amplitudes of different non-neuron cell types were proportional to their estimated proportion in healthy samples (Figure 4b), suggesting that they were likely to be the passive consequence of neuron proportion decrease.
To further investigate the contribution of composition differences in the autistic brain transcriptome alteration, we dissected the bulk gene expressions of all the PFC samples, including the 34 ASD cases and 30 healthy controls within the same sample age range, into the composition-dependent and composition-independent components, based on the DR-based composition pattern. Natural spline based ANCOVA (Methods) was applied to each of the two components to identify genes with DE between ASD cases and healthy controls. Among the 14,032 detected genes, 74% of them (10,358) showed DE in their composition-dependent components (Benjamini-Hochberg corrected FDR<20%), while only no more than 5% of the genes (644) showed DE in composition-independent components under the same criteria (Supplementary Table  S4). The two numbers based on CIBERSORT-based composition pattern were 8,030 and 275, with majority overlapping with the DR-based ones (Supplementary Figure  S4). Both the two sets of genes were significantly overlapping with the 1,775 DE genes identified based on bulk expression (referred as bulk DE genes, Figure 5a, Fisher's exact test, P<0.0001), suggesting that both of them were important contributors in transcriptome alteration in ASD brains.
Although majority of the transcriptome alterations in ASD brains were related to the cellular composition changes, the alterations in the composition-independent components of brain transcriptome were of particular interest. They might represent the active regulatory disruptions in particular brain cells. The functional enrichments (DAVID 31 ) suggested that they were enriched for membrane and synaptic proteins ( Figure 4b, Supplementary Table S5). Consistently, those genes showed significant expression enrichment in adult and fetal quiescent neurons ( Figure 4b). More interestingly, the genes with significant alterations in their composition-independent components, but not those with significant alterations only in their composition-dependent components, contained a significant excess of genes with genetic variants associated with ASD, based on analysis of genes extracted from SFARI AutDB 38 (Fisher's exact test, odds ratio=1.52, P=0.009, Figure 4b), as well as genes from AutismKB database 39 (Fisher's exact test, odds ratio=1.37, P=0.008, Figure 4b). Furthermore, those genes were also consistently identified as differentially expressed genes in autistic brains in previous studies 12, 13 (Fisher's exact test, P<0.005, Figure 4b). It is also worth to note, that 40 microRNA (miRNA) families showed enrichment of predicted binding at the 3'UTRs of those genes (TargetScan v7.1 40 ; hypergeometric test, P<0.05), which was significantly more than by chance (permutation test, P<0.01, FDR=0.65). Some of the enriched miRNAs, including has-miR-363-3p which was also reported to be dysregulated in Alzheimer's disease 41 , have been reported as differentially expressed in autistic brains 42 . All the above observations suggested that the disruptions of neuronal transcriptional regulations in autistic brains, which were likely to be partially caused by the dysregulations of certain neuronal miRNAs, may result in pathological changes in synapses and largely contribute to ASD pathology.

Discussions
In this study, we attempted to investigate the role of changes related or independent to cell type composition in human brain development and pathogenesis of ASD using computational approaches. With the gene expression of eight main cell types in human brain estimated using the human brain single cell RNA-seq data 21 , we applied a two-step transcriptome deconvolution procedure to the human brain postnatal age series RNA-seq data 11 and the ASD brain RNA-seq data 14 . Based on the estimated composition patterns, gene expression was dissected into the composition-dependent and independent component, allowing the further investigation on each of them.
The transcriptome deconvolution in healthy samples resulted in temporal composition patterns consistent with the prior studies. On the other hand, we also notice that our estimated proportion of total neurons in the adult samples reached around 60-70%, which is much higher than the estimation based on cell counting 43, 44 or DNA methylation deconvolution 45 . This may be explained by the fact, that the estimated composition of transcriptome deconvolution actually reflects the proportion of RNAs contributed by each cell type instead of the proportion of cell numbers. In fact, the RNA content in neurons is reported to be close-to-two-fold as much as in glia cells 46 .
Based on the variance analysis, we regard the composition change to be the main source of the age-related expression change in human brain. This is a consistent conclusion drawn by previous studies based on DNA methylation 45 . The proportion of variance explained by cell types, on the other hand, is smaller: only about 30% for the age-related expressed genes, while it is around 50% for the DNA methylation. Although the influence of sampling and measurement noise could not be ruled out, such observation may imply the additional regulatory signal that show age-related manner and is independent from the DNA methylation, which might be an interesting focus for further study.
Interestingly, genes with age-related changes in the composition-dependent components show high expression level in the infant stages, i.e. < 2 years old, with significantly enriched expression in the three types of neurons. This suggests that although some other cell types such as oligodendrocytes and OPCs presented mixing proportion changes in the age-related manner, most of the composition-dependent changes are due to the transition from fetal neurons to adult neurons, i.e. the neuron maturation process. It therefore implies that the neuron maturation is the primary factor creating the age-related expression, especially during the infant stage.
On the other hand, the age-related changes to the composition-independent component is appealing. Unlike the genes with age-related composition-dependent changes which mostly express highly in the early postnatal development, genes with age-related composition-independent changes tend to have higher expression in adults. What's more, these composition-independent changes are significantly related to either synapse in neurons (G1), or extra-cellular regions and signal peptides in astrocytes and endothelial cells (G3), both of which are relevant to cell-to-cell communications. This is unlikely to be a coincidence. As an information processing system, such communications are critically important for human brains 47, 48 , and our results suggest that the complexity growth of the communication system not only depends on the increased number of computational units, i.e. neurons, but also greatly relies on the enhanced inter-cellular communications which are independent from the cell type composition changes across lifespan. Such communications may not limit to the synaptic connections between neurons, but also include the neuron-glia and glia-glia communications which are also critical to the neuronal network functions 49 .
In autistic brains, the contradictory statements of unchanged density of both neuronal and glial cell pools 50 and even more neuron numbers 51 has been reported. Therefore, it is a surprising observation that the estimated neuron proportions in autistic brains are significantly smaller than that in healthy brains. While the influence of unknown technical issue cannot be ruled out, this result implies the possibility that the general transcriptional activity of neurons is inhibited in ASD. In such case, the total mRNA content in neurons reduces, therefore resulting in smaller contribution of neuronal mRNAs to the bulk mRNA pool. This assumption is supported by the depletion of neuronal marker MAP2 in autistic brains 50 . However, further experimental measurement is also awaited for validation.
Considering that neuroinflammation is thought to play important roles in ASD pathology 52, 53 , it is also unexpected to find no evidence suggesting general activation of microglia, the macrophage cells participating neuroinflammation in the central nervous system. Although the increase of microglia proportion is observed, no significant excess of increase is observed when comparing to other glial cell types after normalizing to their relative proportion in healthy brains. Nevertheless, our results did not rule out the importance of neuroinflammation in ASD, considering that the strongest activation of microglia in ASD patients happens in cerebellum and white matters 53, 54 . Although higher microglia density was also reported in ASD dorsolateral PFC 55 , they may not be fully activated, making our analysis lack of power to capture such signal. More investigation is required to comprehensively study the role of PFC microglia in ASD brains.
While our results suggest that majority of the transcriptome alteration detected in PFC of ASD patients are passive alteration due to the changes of relative mRNA contribution from cell types, a group of genes attracts most of the focus. Those 644 genes, enriching for membrane protein genes and synaptic functions in neurons, represent the transcriptome alterations which are at least partially independent to the composition changes. This is an interesting observation, suggesting the pathological changes happened to the regulatory network in neurons regulating synaptic functions in ASD. This linkage is further supported by their enriched genetic association with ASD. It is also concordant with previous studies, which report the abnormalities of synaptic functions as well as neuronal dendrites and spines in autistic humans and autistic mouse models 56, 57 . More comprehensive investigations are awaited for the roles of synaptic disruptions in ASD, which may lead to effective pharmacological treatment in the future.
Although our framework paves a way for more exhaustive analyses regarding the contribution of cell type composition to the transcriptome changes, we are well aware of the limitations of our method. For instance, our analysis failed to squarely pinpoint the one or several cell types with the greatest influence on either of the two components, though cell-type enrichment analysis was capable to give a glance. Our analysis greatly relied on the accurate transcriptome measurement of all or at least close to all of the cell types in the bulk tissue, which limits its applications. The experimental validation is relatively difficult. Nevertheless, the fast-developing single-cell technology, especially the single-cell RNA-seq, should be able to comprehensively, simultaneously and separately quantify changes happened in cell type composition, cell type RNA content and cell type molecular profiles. Indeed, the single cell RNA-seq technology has been widely used to describe cell type compositions in different organs including brains 21, 22, 58 , as well as to characterize the molecular processes driving neurogenesis and somatic reprogramming to neurons 59, 60 . The further application of this technology to characterize human brain development, ageing, and pathology of neural disorders including ASD would be sufficient to verify our observation and of great value to comprehensively understand those important processes.
If nothing else, we hope that our attempt to dissect expression into composition-dependent and composition-independent components will inspire further studies to elucidate transcriptome changes in a more comprehensive manner.

Data
The human brain single cell RNA-seq data was retrieved from SRA (SRP057196). The human postnatal age-series brain RNA-seq data with 40 samples was retrieved from GEO (GSE51264). All the RNA-seq reads were mapped to the human genome hg38 with STAR 2.3.0e using the default parameters. The number of reads covering the exonic regions of each protein-coding gene annotated in GENCODE v21 was counted and normalized using the R package DESeq2 for each data set separately. The autistic brain RNA-seq data retrieved from GEO (GSE59288), were processed in the same way. The pre-calculated RPKM of the fetal human brain samples were downloaded from BrainSpan (http://www.brainspan.org/static/download.html).

Deconvolution for cell-type composition
Three different strategies were used for deconvolution for cell-type composition. The first method, quadratic programming (QP) based deconvolution, was to model the gene expression of each cell type signature genes in the bulk tissue sample as a linear combination of its expression in each cell type according to the cell type mixing proportion. Thus, the deconvolution problem for each bulk tissue sample was represented as a constrained linear least-square problem, which was: Here, f was the vector of cell type mixing proportion, and C was the matrix of gene expression of the cell type signature genes in each cell type, while x was the known expression level of the cell type signature genes. This model was widely used in the deconvolution problem 23, 25 , and can be solved using quadratic programming 26 .
The second method, CIBERSORT 24 , which was developed by Newman et al. in 2015, was also used. Instead of the solving the least square problem, CIBERSORT estimates the relative cell type composition based on nu-support vector regression, minimizing a linear ε -insensitive loss function with L2 regularization.
The third method, namely diffusion ratio (DR) based deconvolution, was based on the simple assumption that the expression of a cell type signature gene in the bulk tissue can be seen as its expression in the cell type scaled by the cell type's mixing proportion. It is similar as the recent algorithm MCP-counter 27 but was implemented in a simpler way. Under this assumption, the proportion of cell type i (represented as f i ) was simply estimated as: Here, g represented each cell type signature gene of cell type i. After calculating f i for all the cell types, f was normalized so that ∑ ݂ ൌ 1 .

Deconvolution for cell-type expression profile re-estimation
The similar model as the QP-based deconvolution for cell-type composition described above was used for the second deconvolution to re-estimate cell-type expression profile, which was used for estimation of deconvolution performance and variance analysis. Similarly, the deconvolution problem for each gene can be seen as a constrained linear least-square problem, that is: Here, c was the vector of re-estimated gene expression level in different cell types, and x was the vector of gene expression across the bulk tissue samples. F was the composition matrix with each row representing the estimated mixing proportion in the corresponding bulk tissue sample. This problem was also solved by using quadratic programming as described above.

Performance measurement of transcriptome deconvolution
Two different types of deconvolution performance measurements were used to compare different deconvolution methods.
The first type of performance measurement, namely similarity measurement, aimed to measure the general similarity between the real bulk tissue gene expression and the predicted gene expression based on the cell type gene expression level and the cell type mixing proportion. It is the commonly used measurement of deconvolution performance. In this case, two measurements were applied: Pearson's correlation coefficient (PCC), and root-mean-square deviation (RMSD) which was calculated as the root mean squares of log10-transformed fold changes between the observed FPKMs and the predicted FPKMs (products of estimated cell-type proportions and re-estimated cell type expressions) across all genes.
The second type of performance estimation, namely discrimination measurement, was established based on the following assumption: the accurate deconvolution can result in re-estimated cell type expression estimations which are not only similar to the measured cell type expressions, but also able to recover the cell type identities discriminating one cell type from the others. Using PCC as the proxy of similarity between two cell type expression profiles, the discriminating recovery score (DRS) of cell type i was defined as: Here, e i was the vector of observed log10-transformed FPKM of cell type i, and ො was the vector of re-estimated log10-transformed FPKM of cell type i. The mean DRS across different cell types was then used as the proxy of the overall DRS.

Dissection of composition-dependent and composition-independent expression components
With X denoted as the N × M gene expression matrix with N genes and M samples, and the d × M estimated composition matrix F by any of the mentioned transcriptome deconvolution method with d cell types, the re-estimated cell type expression matrix C with N rows and d columns was obtained as described above. The composition-dependent component was then represented as the matrix X d = CM, while the composition-independent component was represented as X i = X -X d .

ANCOVA based on natural spline with variable degree of freedom
To identify age-related changes, for each gene, with its observed expression level, composition-dependent component, or composition-independent component as the response variable, the age test, an ANCOVA employing the F-test was used to compare the null model: a linear model only with intercept, to a series of alternative models: the natural spline with degree of freedom from two to eight in response to square root transformed sample ages (sqrt-age). The best alternative model was chosen by applying the adjusted r 2 criterion 8 . Genes with BH corrected P<0.05 were considered as genes with its expression, composition-dependent or composition-independent component changed in the age-related manner. For calculating the proportion of variance explained by age, the natural spline model with degree of freedom equaling to eight in response to sqrt-age was used.
The similar ANCOVA framework was also applied to identify changes between autistic samples and healthy samples, by considering the sample ages. The null model was defined as the unified natural spline model with no discrimination of autistic and healthy samples, while the alternative model with varied coefficient in the two groups was compared by employing the F test. Genes with BH corrected P < 0.2 were considered as genes with composition-dependent or composition-independent component changed between the two groups of samples. Similarly, sqrt-age was used as the independent variable.

Cell type expression enrichment analysis
The cell type enrichment analysis was performed based on the average RPKM of each gene in the eight major cell types in human brains including fetal replicating neurons, fetal quiescent neurons, adult neurons, astrocytes, oligodendrocytes, endothelial cells, microglia and oligodendrocyte precursor cells (OPCs), based on the human brain single cell RNA-seq data 21 . The expression specificity of one gene in one cell type was represented by the log10-transformed fold change between its expression level in the cell type and its average expression level in the other cell types. For each cell type, a one-sided Wilcoxon's rank test was applied to compare the expression specificities of genes in the gene list, to the other genes that were detected in both data. 1 . H  u  f  f  m  a  n  K  .  T  h  e  d  e  v  e  l  o  p  i  n  g  ,  a  g  i  n  g  n  e  o  c  o  r  t  e  x  :  h  o  w  g  e  n  e  t  i  c  s  a  n  d  e  p  i  g  e  n  e  t  i  c  s  i  n  f  l  u  e  n  c  e  e  a  r  l  y  d  e  v  e  l  o  p  m  e  n  t  a  l  p  a  t  t  e  r  n  i  n  g  a  n  d  a  g  e  -r  e  l  a  t  e v  a  s  c  u  l  a  r  c  e  l  l  s  o  f  t  h  e  c  e  r  e  b  r  a  l  c  o  r  t  e  x  .   J  N  e  u  r  o  s  c  i   3  4  ,  1  1  9  2  9  -1  1  9  4  7  (  2  0  1  4  )  .   1  6  .  A  l  l  e  n  N  J  ,  B  a  r  r  e  s  B  A  .  N  e  u  r  o  s  c  i  e  n  c  e  :  G  l  i  a  -m  o  r  e  t  h  a  n  j  u  s  t  b  r  a  i  n  g  l  u  e  .   N  a  t  u  r  e   4  5  7  ,  6  7  5  -6  7  7  (  2  0  0  9 ) .  e  l  v  i  g  D  P  ,  P  a  k  k  e  n  b  e  r  g  H  ,  S  t  a  r  k  A  K  ,  P  a  k  k  e  n  b  e  r  g  B  .  N  e  o  c  o  r  t  i  c  a  l  g  l  i  a  l  c  e  l  l  n  u  m  b  e  r  s  i  n  h  u  m  a  n  b  r  a  i  n  s  .   N  e  u  r  o  b  i  o  l  A  g  i  n  g   2  9  ,  1  7  5  4  -1  7  6  2  (  2  0  0