Marcel Leist

Philipp Bastian Zimmer

Balthasar Kuegler

Birte Baudis

Andreas Genewsky

Vivek Tanavde

Winston Koh

Betty Tan

Tanja Waldmann

Suzanne Kadereit

Coordinated waves of gene expression during neuronal differentiation of embryonic stem cells as basis for novel approaches to developmental neurotoxicity testing

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Pluripotent stem cells are able to form any cell type in an organism, including all different types of neurons. Ultimately, the entire complexity of the mammalian central nervous system (CNS) is generated during ontogenesis from a few single cells. This intricate process involves proliferation and cell fate determination events as well as migration, synaptogenesis, apoptosis and myelination. Many of these steps, in particular neuronal generation and differentiation, can be recapitulated by embryonic stem cells (ESC) under appropriate culture conditions [START_REF] Conti | Neural stem cell systems: physiological players or in vitro entities?[END_REF][START_REF] Gaspard | An intrinsic mechanism of corticogenesis from embryonic stem cells[END_REF][START_REF] Kuegler | Markers of murine embryonic and neural stem cells, neurons and astrocytes: reference points for developmental neurotoxicity testing[END_REF][START_REF] Barberi | Neural subtype specification of fertilization and nuclear transfer embryonic stem cells and application in parkinsonian mice[END_REF][START_REF] Götz | The cell biology of neurogenesis[END_REF][START_REF] Abranches | Neural differentiation of embryonic stem cells in vitro: a road map to neurogenesis in the embryo[END_REF] .

ESC-based studies of neurodevelopment allow investigations and interventions not

easily possible in vivo. However, known differentiation protocols differ in their suitability for such mechanistic or even toxicological studies. They were originally developed to generate large cell numbers e.g. for cell substitution therapies, and often focus on one defined neuronal subtype [START_REF] Barberi | Neural subtype specification of fertilization and nuclear transfer embryonic stem cells and application in parkinsonian mice[END_REF][START_REF] Lee | Efficient generation of midbrain and hindbrain neurons from mouse embryonic stem cells[END_REF] . For instance, older protocols involve a step of embryoid body (EB) formation before neural induction is triggered [START_REF] Strübing | Differentiation of pluripotent embryonic stem cells into the neuronal lineage in vitro gives rise to mature inhibitory and excitatory neurons[END_REF] . EBs contain cells of all three germ layers, and neurogenesis is often enhanced by addition of the morphogen and patterning factor retinoic acid (RA) [START_REF] Lee | Efficient generation of midbrain and hindbrain neurons from mouse embryonic stem cells[END_REF] . Frequently, only a small number of the initially-present ESC forms neurons, and the observation of individual cells is hardly possible. Other protocols use co-cultures with stromal cell lines like MS5 [START_REF] Barberi | Neural subtype specification of fertilization and nuclear transfer embryonic stem cells and application in parkinsonian mice[END_REF] to differentiate ESC towards neurons, and would therefore introduce additional complexity into models for developmental neurotoxicity (DNT). A recently developed monolayer differentiation protocol allows monitoring of the differentiation procedure and of possible effects of different chemicals during the whole period of differentiation on a single cell level. Moreover, the use of fully defined media components allows for stringent system standardization [START_REF] Gaspard | An intrinsic mechanism of corticogenesis from embryonic stem cells[END_REF][START_REF] Abranches | Neural differentiation of embryonic stem cells in vitro: a road map to neurogenesis in the embryo[END_REF][START_REF] Ying | Defined conditions for neural commitment and differentiation[END_REF] .

DNT is the form of toxicity least examined and hardest to trace, as it is not necessarily related to cell loss. Behavioural pathology in the absence of cell loss is also known from disease models, e.g. for Huntington's disease [START_REF] Hansson | Transgenic mice expressing a Huntington's disease mutation are resistant to quinolinic acid-induced striatal excitotoxicity[END_REF] or schizophrenia [START_REF] Penschuck | Decrease in parvalbumin-expressing neurons in the hippocampus and increased phencyclidine-induced locomotor activity in the rat methylazoxymethanol (MAM) model of schizophrenia[END_REF] . Toxicants, such as mercury, lead or polychlorinated biphenyls may trigger behavioral or cognitive deficits without histophathological hallmarks [START_REF] Grandjean | Developmental neurotoxicity of industrial chemicals[END_REF] . As the DNT endpoints are particularly difficult to test, less than 0.1% of frequently used industrial chemicals have been examined, and for the few known toxicants the mechanism of action is still elusive (reviewed in [START_REF] Grandjean | Developmental neurotoxicity of industrial chemicals[END_REF][START_REF] Makris | A retrospective performance assessment of the developmental neurotoxicity study in support of OECD test guideline 426[END_REF][START_REF] Bal-Price | In vitro developmental neurotoxicity (DNT) testing: Relevant models and endpoints[END_REF] ). It has been suggested that cellular physiology (Ca 2+ handling) may be affected during the period of exposure [START_REF] Rossi | Modifications of Ca2+ signaling by inorganic mercury in PC12 cells[END_REF] . This may eventually lead to changes in differentiation and patterning in the CNS, which is the basis for long term effects that are observed after the exposure to toxicants has ceased.

CNS development is assumed to be orchestrated by waves of gene expression [START_REF] Wei | Global gene expression patterns during neural differentiation of P19 embryonic carcinoma cells[END_REF][START_REF] Aiba | Defining a developmental path to neural fate by global expression profiling of mouse embryonic stem cells and adult neural stem/progenitor cells[END_REF] that determine different intermediate cell phenotypes and form the basis for subsequent steps. Some periods may be more sensitive to certain toxicants than others. Epidemiological proof for such "windows of sensitivity" in organ development with long term consequences for the organism comes from thalidomide exposure in man [START_REF] Kuegler | Markers of murine embryonic and neural stem cells, neurons and astrocytes: reference points for developmental neurotoxicity testing[END_REF] and various animal models [START_REF] Jongen-Relo | The prenatal methylazoxymethanol acetate treatment: a neurodevelopmental animal model for schizophrenia?[END_REF] .

The use of ESC may allow for new approaches to understand mechanisms of DNT and to evaluate the safety of chemicals. However, current test systems based on the differentiation of stem cells to either cardiomyocytes [START_REF] Marx-Stoelting | A review of the implementation of the embryonic stem cell test (EST). The report and recommendations of an ECVAM/ReProTect Workshop[END_REF] or neural cells [START_REF] Bal-Price | In vitro developmental neurotoxicity (DNT) testing: Relevant models and endpoints[END_REF] neither yield mechanistic info, nor do they account for the complexity of CNS development, i.e. the establishment of a balance between multiple neuronal cell types [START_REF] Kuegler | Markers of murine embryonic and neural stem cells, neurons and astrocytes: reference points for developmental neurotoxicity testing[END_REF][START_REF] Rao | Developmental Neurobiology[END_REF] . The "toxicology for the 21 st century" initiative [START_REF] Collins | TOXICOLOGY: Transforming Environmental Health Protection[END_REF][START_REF] Leist | The dawning of a new age of toxicology[END_REF] suggests the identification of pathways, and proposes a biologically-based explanation of toxicant effects, as opposed to the current black-box test systems. In the case of ESC-based models of DNT, this requires a detailed understanding of the developmental process leading to multiple different cell types.

Detailed knowledge on the waves of gene induction controlling different developmental steps would be an essential prerequisite. However, CNS development is proceeding at different paces. For instance, the anterior and posterior part of the neural tube follow different kinetics, and some regions of the CNS continue neurogenesis, while in other regions cells have already reached fully postmitotic stages of their cells [START_REF] Rao | Developmental Neurobiology[END_REF] . At present, it is not clear to which extent this is replicated in developing ESC cultures, but it is known that waves of neuronal development of different amplitudes may be detected in vitro [START_REF] Wei | Global gene expression patterns during neural differentiation of P19 embryonic carcinoma cells[END_REF][START_REF] Aiba | Defining a developmental path to neural fate by global expression profiling of mouse embryonic stem cells and adult neural stem/progenitor cells[END_REF] . These would need to be deconvoluted and characterized in depth in order to provide markers that comprehensively describe the neurodevelopmental process.

Our study was undertaken to analyze the wave-like expression pattern of mESC neurodevelopment as a basis for the definition of test windows and markers. This knowledge should help to identify non-cytotoxic, but neuroteratogenic compounds able to shift neuronal composition or phenotypes. Finally, the markers should distinguish multiple cell types and differentiation stages, and be able to indicate subpopulations of adherent cells that are inhomogeneously distributed.

Results

Monolayer differentiation of mESC to neurons

CGR8 were chosen as a widely available murine ESC line suitable for feeder-free culture maintenance and with established potential to develop along the neuroectodermal and neuronal lineage [START_REF] Suter | A Sox1 to Pax6 switch drives neuroectoderm to radial glia progression during differentiation of mouse embryonic stem cells[END_REF][START_REF] Conti | Nicheindependent symmetrical self-renewal of a mammalian tissue stem cell[END_REF] . Differentiation was performed under defined serum-free and adherent conditions according to a published protocol [START_REF] Ying | Defined conditions for neural commitment and differentiation[END_REF] . The method was optimised for developmental neurotoxicity testing requirements, which demand high reproducibility between individual experiments. The efficiency of the final differentiation was confirmed by immunocytochemistry on day-of-differentiation 20 (DoD20). The majority of cells was positive for the pan-neuronal markers Tuj1 (neuronal form of betaIII-tubulin) and NeuN (neuronal-specific postmitotic nuclear antigen, encoded by fox3) [START_REF] Kim | Identification of Neuronal Nuclei (NeuN) as Fox-3, a New Member of the Fox-1 Gene Family of Splicing Factors[END_REF] . Many cells also expressed synapse associated markers such as SV2 or PSD95, respectively (Fig. 1A). As a more quantitative overall measure for the robustness of the differentiation protocol to generate neurons from mESC, we chose mRNA expression, which we followed over time in five independent differentiation experiments, performed at different times, with different CGR8 cell batches, and by different operators. The kinetics of loss of the stemness marker oct4, transitional upregulation of the neural stem cell marker nestin and the delayed induction of the mature neuronal marker synaptophysin or the glial marker gfap were highly reproducible across different experiments (Fig. 1B). For functional characterization, neurons from at least three independent differentiations were tested for electrophysiological activity. Differentiation to mature, electrophysiologically-active neurons was shown by the presence of voltagedependent Na + and K + and Ca 2+ channels in individual patch-clamped neurons (Fig. 2A-C, Fig. S1). Further experiments also identified spontaneous neuronal electrical activity (Fig. 2D). Action potentials could be evoked by depolarization of individual neurons (Fig. S1), and currents were also evoked by exposure to N-methyl-D-aspartate (NMDA) or kainic acid. The latter were blocked by the respective selective antagonists of NMDA and non-NMDA glutamate receptors (Fig. 2E). Alternatively, picrotoxin-sensitive currents were evoked by GABA (Fig. S1). Accordingly, the chosen neuronal differentiation protocol provides a solid basis for reproducible generation of bona fide neurons.

Transcription-based endpoints to identify disturbed neuronal differentiation

We next investigated whether we could detect subtle perturbations of the differentiation process below the cytotoxicity threshold. We tested whether mRNA-based readouts would fulfil the requirements of giving such information, e.g. altered neuronal patterning or shifts in subpopulations. Parallel mESC cultures were differentiated for 7, 15, and 20 days, and mRNA was prepared for quantitative RT-PCR analysis. These cells were treated during two different time windows (DoD1-7, DoD8-15) with two neuro-teratogens (Fig. 3A). With the concentrations used here cell death was not detectable (data not shown) and cells looked viable and were morphologically indistinguishable from untreated cells (Fig. 3B). We used the morphogen retinoic acid (RA) as a known in vivo and in vitro reproductive toxicant and cyclopamine for its ability to alter sonic hedgehog (Shh) signaling resulting in the disruption of patterning gradients responsible for floor plate and ventral neurons [START_REF] Gaspard | An intrinsic mechanism of corticogenesis from embryonic stem cells[END_REF][START_REF] Rao | Developmental Neurobiology[END_REF] . The analysis of neural patterning-associated marker genes in the chemical-treated cultures showed that their differentiation had been severely affected compared to untreated cultures. As expected from the literature [START_REF] Irioka | Distinct effects of caudalizing factors on regional specification of embryonic stem cell-derived neural precursors[END_REF] , RA induced accelerated neuronal differentiation indicated by increased synaptophysin (marker of mature neurons) expression (Fig 3A,a,b) whereas cyclopamine reduced the expression of markers typical for more ventrally-located neurons like Shh, Nkx2.1 and Dlx1 (Fig 3A,a) but not overall neuronal differentiation (Fig. 3A). These data indicate that well chosen marker genes can indicate subtle shifts in differentiation patterns not visible morphologically. This was further analyzed in different variations of the incubation scheme with a set of stage-specific markers (Fig. 3A, b-d). For instance cultures were exposed to cyclopamine from DoD1-7 (Fig. 3A,b) and immediately analyzed thereafter. Treatment did not affect the overall formation of NPC (as indicated by unchanged nestin expression), but again the reduced Shh expression suggested a reduced ventral development. When cells were released from cyclopamine after 7 days of treatment and left to differentiate further the compound effect on shh expression was still observable on day 15 (Fig. 3A,c). The known dorsalizing activity of cyclopamine would suggest a shift of neurotransmitter phenotype from GABAergic (Gad2 as marker) to glutamatergic (vglut1 as marker) [START_REF] Gaspard | An intrinsic mechanism of corticogenesis from embryonic stem cells[END_REF] . This was not observed after treatment for the first seven days (Fig 3A,c), but a significant decrease of Gad2 (more ventrally prominent) was observable when the cells were treated during a short transient period between DoD8 and 15 (Fig. 3A,d).

In the case of RA, the acceleration of development (increased synaptophysin expression) was already significant at early stages (Fig 3A,b). In addition we examined its known caudalizing effect [START_REF] Irioka | Distinct effects of caudalizing factors on regional specification of embryonic stem cell-derived neural precursors[END_REF] . We found upregulation of markers usually expressed in caudal parts of the neural tube (hoxa6, hb9), and associated with the development of motor neuron precursors (isl1) (Fig. 3A,c,d).

Finally, we examined whether inhibited (instead of shifted) differentiation by non-cytotoxic chemicals was detectable by RNA markers. We used 3i, a kinase inhibitor mix known for inhibiting differentiation of mESC [START_REF] Ying | The ground state of embryonic stem cell self-renewal[END_REF] . Early exposure (Dod0-7) resulted in cultures with an immature cellular phenotype (data not shown) and retarded neural differentiation indicated by a decreased expression of neural markers hes5, nestin and betaIIItubulin and an increased expression of the stemness factor Oct4 (Fig. S2,a). When cells were treated on DoD1-7, washed thoroughly, and then left to differentiate without interference until DoD15, 3i did not show effects on Oct-4 expression anymore, but neural differentiation was severely delayed as seen by significantly lowered expression of nestin, hes5 and the mature neuronal marker synaptophysin (Fig. S2,c). Treatment of cells with 3i after DoD7 (after neural differentiation had been initiated) did not return them to the stem cell state, but was rather cytotoxic. In summary the usefulness of transcript profiling for detection of patterning disturbances was clearly demonstrated and confirmed.

The low variation of the data across different experiments suggests that transcript profiling is useful for quantitative assessments of disturbed neuronal differentiation. As the effects of chemicals are likely to be different, depending on the time window of exposure and on the differentiation stage assessed, it may be necessary to measure the impact on differentiation at different DoD, and specific markers need to be selected for each stage. In order to identify such potential markers, global, and kinetically-resolved transcription profiling was performed.

Identification of clusters of genes regulated during neuronal differentiation of mESC.

We analyzed changes in the transcriptome over time (DoD 0, 7, 15, 20) based on oligonucleotide microrarrays. By using a few marker genes, we verified that the differentiation kinetics of the cultures used for microarray analysis matched the ones observed during many other well-controlled experiments performed earlier (Fig. 1B). Having established this important correlation, the kinetics of expression of each gene represented on the chip was used as input for an unbiased clustering analysis. This analysis yielded eight groups of regulation profiles (Clusters Ia, Ib, IIa, IIb, IIIa, IIIb, IV and V) (Fig. 4A,S3), besides the group of genes not regulated at all or not yielding any signal in this analysis (genes not shown). Cluster Ia was characterized by rapid (within the first 7 days), and cluster Ib by slow (continuous from DoD0 to DoD20) downregulation. These two clusters obtained by non-hierarchical analysis exemplify the principle of superimposed gene regulation waves with different amplitudes (see also clusters III-V). Clusters IIa and IIb contained genes that were transiently regulated at DoD7 (Fig. 4A) (IIa: transient upregulation, IIb: transient downregulation). Cluster IIIa and IIIb were characterized by a rapid increase of transcripts between day 0 and DoD7 maintained then at high levels (Fig 4A). Cluster IV contained genes, which remained low until DoD7 and then reached high levels on DoD15, with no major changes until DoD20 (Fig 4A). The final cluster V comprised transcripts that were hardly upregulated until DoD15, and reached their maximum on DoD20 (Fig. 4A). Thus, large groups of genes appear to be regulated in a co-ordinated manner in defined, and partially overlapping, waves. Notably, the most-pronounced changes in expression occurred during culture periods not involving any manipulation of the cells, i.e. between DoD0 and DoD7 (cells are plated on day 0 and only replated on day 7), and between DoD15 and DoD20. This has important implications for the choice of relevant treatment windows in studies of developmental neurotoxicity.

To explore this further, the genes falling into the different clusters were subjected to a more detailed analysis. First, we controlled for the efficiency of differentiation, and examined the behaviour of about 40 genes that characterize the initial mESC stage 3 . We found 33 of them differentially regulated in our analysis and examined their distribution across the individual clusters. About 80% of these fell into cluster Ia and Ib (overall downregulation), while the remainder was found in cluster IIb (downregualtion from DoD0-DoD7). Thus, all mESC markers identified on the chip were confirmed to be downregulated upon initiation of the neuronal differentiation (Fig. 4B). Next, we proceeded to ask how neural precursor cell (NPC) and neuronal differentiation (N) markers clustered. About 130 N-markers were found in the different clusters. NPC markers (n = 73) were assembled from the literature 3 , and 63 of them were found in our analysis. Concerning the relative distribution between the clusters, most NPC markers were found in clusters IIa and IIIa/b (both containing genes upregulated early) and about 15% were found in clusters IV and V (Fig. 4B). In contrast to this, most N markers were found in the clusters with increasing gene expression (clusters III-V), while about 20% were found in cluster IIa (transient upregulation on DoD7) (Fig. 4B). These findings strongly suggest that the clusters identified by unbiased bioinformatics methods, correspond to waves of real biological processes describing the differentiation process of mESC to neurons. To explore this working hypothesis further, we continued with a detailed analysis of the biological significance of genes in individual clusters.

Loss of pluripotency is accompanied by progressive changes in transcripts responsible for chromatin organization and DNA/cell cycle functions

Initially, genes in cluster I were analyzed for GO categories significantly overrepresented. For instance, cell cycle-related processes were identified (Fig 5A ,S4). More unexpectedly, also chromatin structure and epigenetic processes seemed to be affected (Fig 5A). We examined this in more detail. Gene lists of relevant processes were assembled both with the help of the GO data base and extensive literature search. The clusters were then queried for the presence of these genes (Fig. 5A,S4). Interestingly, of 14 genes known to be associated with chromosome structure, all were identified in clusters Ia and Ib (overall downregulation). The same applied to genes associated with DNA replication (n = 57), DNA repair (n = 62) and DNA methylation (n = 5). Also, most of the genes coding for histones, histone modifiers, chromatin remodelling and chromatin substructuring were found in cluster I (Fig. 5A,S4).

We also examined whether the transcriptional changes gave clear evidence for global quantitative changes in chromatin structure. However, the observed pattern of regulated genes rather suggested a restructuring of the chromatin without an overall increase of DNA methylations or histone modifications (Fig. S4). Therefore restructuring of pericentric heterochromatin was examined by high resolution confocal microscopy. While the chromatin, distributed relatively homogeneously over the nucleus in mESC, it was organized entirely differently after 20 days of differentiation, with neurons characterized by conspicuously large (1-2 µm scale) dense areas of chromatin (Fig. 5B). In general, many genes of cluster Ia/b may be typical for the biological processes initiated in an actively cycling stem cell developing towards a post-mitotic and highly differentiated cell. As indicated by our data, this involves many different processes running in parallel. Among the identified genes, four (Smarca1, Myst4, Jmjd3 and Hdac11) are known to be neurospecific, and five (Suz12, Ezh2, Bmi1, Cbx2 and Cbx8) are components of the polycomb repressor complexes (PRC), which play an important role in differentiation-related control of gene promoters. These genes could serve as sensitive markers to detect negative effects of compounds on early developmental processes (Fig. S4).

Correlation of neural precursor formation with a strong, transient change of gene expression levels

We examined whether genes of cluster II (Fig. 4A) were specifically linked to the process of neural precursor cell (NPC) formation, an essential step in neurodevelopment. One of the genes to be expected, and actually found, in cluster IIa was nestin(Fig. 1B), an established marker for NPC [START_REF] Kawaguchi | Nestin-EGFP transgenic mice: visualization of the self-renewal and multipotency of CNS stem cells[END_REF] . We confirmed nestin expression by immunocytochemistry, with a majority of cells expressing this cytoskeletal protein on DoD7. Moreover, nestin-positive cells often seemed to be arranged in ring-like structures, reminiscent of so-called rosettes, structures resembling two-dimensional neural tubes [START_REF] Elkabetz | Human ES cell-derived neural rosettes reveal a functionally distinct early neural stem cell stage[END_REF][START_REF] Zhang | In vitro differentiation of transplantable neural precursors from human embryonic stem cells[END_REF] (Fig. 6A). Quantification by flow cytometry analysis showed that about 80% of all cells in the culture became nestin-positive (Fig. 6B).

High synchronization at this stage was suggested by the very distinct and sharp expression profile of genes in cluster IIa (Fig. 6C). Some genes were upregulated several thousand-fold on DoD7, and strongly downregulated again on DoD15 and DoD20. Besides nestin, many other genes typically associated with neuroepithelial precursors, NPC, and neurogenesis were found in cluster IIa (S3). In addition, some genes associated with early, but definitive neuronal development, and with specification and patterning of neuronal subtypes were identified (e.g. Dll1, Hes3). However, cluster IIa also contained genes thought to be involved in multiple cellular and metabolic processes besides neurodevelopment (e.g. Janus kinase 2 (Jak2), forkhead box D4 (Foxd4), B-cell leukemia/lymphoma 2 (Bcl-2), kinesin family member 21A (Kif21a), angiotensin II receptor, type 1a (Agtr1a, also known as AT1), monooxygenase-DBH-like 1 (Moxd1), acetoacetyl-CoA synthetase (Aacs), ADP-ribosylation factor-like 2 binding protein (Arl2bp), stearoyl-Coenzyme A desaturase 2 (Scd2). Therefore, we performed an unbiased and quantitative analysis of the biological significance of cluster IIa, and examined which of the 16000 gene ontology (GO) categories were statistically overrepresented by cluster IIa genes. The result confirmed the efficacy of the chosen differentiation protocol, as the GO "nervous system development" emerged with a p-value < 10 -13 , and only neuronal/neurodevelopmental GOs were identified with the exception of eight (weakly significant) potential ossification genes (Table 1). Thus, genes of cluster IIa are associated with determination of early neuroectodermal/neuronal development and their expression could represent an important endpoint for testing of disturbed proliferation and differentiation during this early and crucial developmental time window.

Markers of regional fate decisions in the CNS

Overall in vivo neurodevelopment is a non linear sequential process, even at this initial stage. Early central nervous system formation is characterized by local differences in timing and patterning [START_REF] Rao | Developmental Neurobiology[END_REF] that may be reflected in some in vitro systems, while others show already a clear regional fate decision 1 . To characterize our differentiation system more thoroughly, we examined region-specific markers. Expression of regional markers was relatively limited in cluster IIa. Amongst the few markers expressed, forebrain (e.g. Foxg1) and hindbrain (e.g. Hoxa2/b2)-related indicator genes were evenly distributed in cluster IIa (Fig. 6D,S5). More patterning markers were found in the clusters containing continuously upregulated genes (clusters IIIa/b and IV+V). Also here, forebrain (e.g. Reln), midbrain (e.g. En1/2), and hindbrain (e.g. Lmx1a or Hoxa1) markers were evenly distributed (S5). Accordingly, our experimental model appears to reflect several parallel lines of in vivo neural specification. In this system with broad developmental potential, the ratios of expression of different patterning markers may provide very sensitive indicators of disturbed neurodevelopment worth an in depth exploration as piloted in figure 3.

Specificity for neuronal induction with respect to glial cells

Most neuronal differentiation systems generate a certain percentage of neurons, but frequently also less-characterized populations of additional neural cells plus non-neural cells.

Characterization of these undesirable cell types within the culture is of particular importance for quantitative assays of disturbed neurodevelopment, where such populations may increase on cost of neurons. Moreover, such cells could confound (enhance or decrease) neurotoxicity of tested compounds. The now available broad transcriptional profile allowed us a more detailed analysis. First, we examined gliogenesis, as Gfap was sharply upregulated on DoD20 (Fig. 1B), and some small GFAP-positive cell areas were reproducibly identified by immunocytochemistry (Fig. 7A). These cells displayed typical astrocytic morphology and were found in 1-2 small islands/cm 2 . An unbiased search for overrepresented GO categories did not result in any hits related to gliogenesis or glial function (Table 2). As alternative approach, we used a list of 25 astrocyte-related genes [START_REF] Kuegler | Markers of murine embryonic and neural stem cells, neurons and astrocytes: reference points for developmental neurotoxicity testing[END_REF] and found 11 of them to be upregulated on DoD20 compared to DoD0, with 4 additional astrocytes-related genes transiently upregulated on DoD7 (Fig. 7B). Most likely, the early upregulation (e.g. vimentin) of apparent astrocytic markers is due to the generation of radial glia-like NPC at DoD7. This cell type, as exemplified by the upregulation of Fabp7 (brain lipid binding protein; in cluster IIIb) [START_REF] Feng | Brain lipid-binding protein (BLBP): a novel signaling system in the developing mammalian CNS[END_REF] or ascl1 (=Mash1) [START_REF] Battiste | Ascl1 defines sequentially generated lineage-restricted neuronal and oligodendrocyte precursor cells in the spinal cord[END_REF] shares many markers with astrocytes [START_REF] Götz | The cell biology of neurogenesis[END_REF] . In contrast to this, the late upregulation of many mature astroglial markers is due to a real generation of astrocytes, as corroborated independently by immunocytochemistry (Fig. 7A) and by qPCR of a larger set of markers and (Fig 1B ,7C). Small increases of this normally minor cell population may affect toxicity testing during the later differentiation phases and requires good controls on the basis of mRNA markers identified here. Any potential role of astrocytes appears to be insignificant for early phases of DNT testing.

We also examined genes known to be expressed in oligodendrocytes and microglia. The latter cell type appeared to be absent. Few oligodendrocyte markers were upregulated, while many functionally important genes were not (Fig. 7B). By immunocytochemistry, we were not able to identify any oligodendrocytes either (not shown). Thus their contribution to the DoD20 cultures appears to be negligible.

Specificity for neuronal induction with respect to other germ layer lineages

Finally, all GO categories significantly overrepresented by the genes of clusters III-V (upregulation on DoD20 vs DoD0) were determined bioinformatically, and searched for evidence of non-neural cell type formation. Individual clusters did not indicate any non-neural cell types while representation of neuronal GOs was highly significant (Table 2). Upon pooled analysis of clusters IV and V (late upregulation), the GOs "blood vessel development" and "muscle organ development" emerged as significant. Thus, it is possible that a subpopulation of cells present on DoD20 may have a smooth muscle phenotype or partial features of such cells.

Waves of clustered genes related to neuronal induction

The classic characterization of neuronal differentiations is heavily based on immunocytochemistry (Fig. 1A). This approach allows only for detection of a few antigens, and quantitative analysis has been difficult in cultures with heterogeneous cell distribution.

Here, we used instead transcriptional profiling. For characterization of the cultures, we used groups of genes identified by non-biased clustering as basis to bioinformatically extract the main GO categories represented by them (Table 2). In addition, a complementary approach based on literature and expert judgement was chosen to select interesting groups of genes (Fig. 8). Novel information emerged from this combination of strategies. Most importantly, we found that the differentiation did not proceed as a linear sequence of synchronized sequential steps, although clear and distinct waves of genes had been identified by cluster analysis. The process rather involved strongly overlapping processes with one underlying large wave (cluster IIIa/b) superseded by shorter waves (cluster IV and V). For instance, generation of neurons and axogenesis/growth cone formation seemed to be ongoing in the entire period from DoD7 to DoD20 as indicated by groups of neuroreceptors and growth cone/axon guidance-related genes in cluster III (Fig. 8A,B). A larger group of genes associated with synaptic vesicles or the transmission of nerve impulse only appeared later (cluster IV/V) (Fig 8A). In the latest phase, hardly any new genes associated with growth cones and neurodevelopment were induced. Instead, genes associated with responses to stress and hormonal stimuli were strongly up-regulated (Table 2). Interestingly, one frequently neglected group of genes also emerged significantly: the "regulation of extracellular matrix components" (Fig 8C). Such genes were strongly and distinctly up-regulated in late differentiation stages, similar to the more expected transporters, genes associated with neuronal projection and post-synaptic membranes. A further interesting group of genes was highly upregulated during the last step of differentiation and involved those known to be "associated with hereditary neurodegenerative diseases". Reanalysis by quantitative PCR confirmed this finding and yielded regulation factors (DoD20 vs DoD0, n = 2 differentiations) of: 92-fold and 16-fold for the Alzheimer's disease associated genes App and Mapt, 273-fold for the schizophrenia-associated gene Nrnx1, 91-fold for the prion protein Prnp, and 19-fold/ 56-fold for the Parkinson's disease-related genes Pink1/ Snca. We wondered whether toxicants may affect this very late phase of neuronal differentiation identified by our transcriptional profiling approach. We chose lead as well-established toxicant with a purported role in the developmental origins of neurodegenerative diseases (see Fig. S6). The cells were treated with a non cytotoxic concentration (assessed by resazurin reduction and LDH release, data not shown) of lead (1 µM) only during the last phase of differentiation (DoD14-DoD20). The transcript levels of two neuronal markers and the set of disease associated genes identified above were used to examine differences in differentiation.

The mRNA expression of Tubb3 and Syp was significantly reduced. Moreover lead exposure had also a dampening effect on the expression of App, Mapt, Nrnx1 and Prnp (Fig. 8D).

Pink1 and Snca were not affected. Also their relative increase with respect to the pan-neuronal marker synaptophysin was not significant. The knowledge on markers of such different processes together with that of the expected timing of their expression provides an ideal toolkit for fine-mapping of subtle developmental disturbances.

DISCUSSION

We have here demonstrated a concept of overlapping waves of gene regulation and suggested its use to define protocols, test windows and endpoints for developmental neurotoxicity testing. Our findings should be helpful to close a gap between two highly developed, but isolated disciplines: experimental developmental neurobiology and toxicology. The former has been highly successful in defining the functional importance, regional expression and cell type/stage association of important genes. The latter has an urgent need for robust and sensitive marker genes to identify disturbances of development. As functionally important genes are not always suitable markers, and good markers sometimes lack a functional role, it was important to examine the complex system of mESC differentiation to neurons to define patterns of gene expression changes useful as toxicological endpoints. We confirmed here earlier findings 2,6,9,27 that the chosen differentiation protocol yields highly reproducible results and supports the development of multiple neuronal subtypes, while maintaining a high specificity for neurons as opposed to other cells. In this system we showed that subtle changes in the speed of differentiation, or in dorso-ventral or anterior-posterior patterning due to toxicants can be detected by using the right choice of mRNA markers. Such changes may be considered in vitro correlates of known teratogenic effects of the chosen compounds. For instance, cyclopamine is known to cause dramatic patterning disturbances in a very defined period of brain development, and a defined brain region (holoprosencephaly); and retinoic acid is well-known to cause shifts in the anterior-posterior axis organisation favouring the more posterior parts, as found here by transcript markers. Lead is known to subtly affect multiple neuronal types, which is in agreement with the broad pattern of disturbances found here for exposure to lead during late neuronal differentiation (see Fig. S6 for references).

However, the data also suggest some warning on the limitations of in vitro -in vivo correlations. Although the cyclopamine data in this study suggest some type of disturbance in patterning, they would not indicate a problem in the separation of the forebrain hemispheres, as observed in animal studies. Similarly, the in vitro test system would not be able to predict the non-neuronal teratogenic effects of retinoic acid. Thus, observations from stem cell systems will have a major value for raising alerts on certain compounds and pinpointing potential mechanisms, while complementary data from other systems may be required to predict specific effects on humans.

Transcriptional profiling studies frequently rely predominantly on stringent bioinformatic analysis. Although this approach appears to be unbiased, it strongly suffers from the weakness and errors of data bases and algorithms. An additional problem is the visualization of the large amount of data in a form that generates meaningful knowledge instead of long tables or unreadable heat maps. In more extensive studies using hundreds to thousands of chips, the generation of data bases with user-friendly interfaces for secondary individual analysis can be a solution e.g. FunGeneES for mESC differentiation [START_REF] Schulz | The FunGenES database: a genomics resource for mouse embryonic stem cell differentiation[END_REF] . This does, however, not solve the inherent bioinformatic problems. For example, assignment of genes to GO categories is not always perfect. For instance, the GO for gliogenesis contains ubiquitous signalling and metabolic molecules (such as Igf1, Sod1) as well as highly specific transcription factors. On the other hand, typical astrocyte markers such as Gfap and glutamine synthetase are not members of this GO. It is evident that bioinformatic analyses to identify markers are problematic on such a basis. First, certain important markers may not be present, second genes may be identified as markers, although they do not qualify biologically and third, the equal weight given to ubiquitous vs. specific genes in statistical analysis results in biological skewing. With these considerations in mind, we chose to combine bioinformatic analysis with classical knowledge-based approaches. The latter were used to control and extend findings based on non-biased approaches, but were also used independently, for the definition and compilation of stem cell and glial markers, and for the analysis of biological functions related to chromatin and neuronal function (Fig 5 and8). During this procedure, the entire hit list of several thousand genes was manually screened, sorted and annotated. For the former approach, a large consortium of experts was consulted, and results were compiled in an open access review format [START_REF] Kuegler | Markers of murine embryonic and neural stem cells, neurons and astrocytes: reference points for developmental neurotoxicity testing[END_REF] . We strongly advocate such combined approaches for toxicological systems biology, which is at present strongly driven by computational methods [START_REF] Leist | The dawning of a new age of toxicology[END_REF][START_REF] Hartung | Food for thought ... on in silico methods in toxicology[END_REF] .

We used different types of endpoints as framework for the transcriptional profiling.

Electrophysiology was chosen as functional endpoint to verify that culture conditions were optimal. The studies show that cells with a distinct neuronal morphology indeed expressed the major types of channels and were electrically-or chemically excitable and able to generate action potentials (Fig 2 ,S1). Although we included replicates from various differentiations, these studies have a rather qualitative character, as the cells that were patched were chosen based on certain morphological criteria and may not represent the entire culture. However, our results fully corroborate earlier findings that functional neurons can be generated from mESC [START_REF] Lee | Efficient generation of midbrain and hindbrain neurons from mouse embryonic stem cells[END_REF][START_REF] Strübing | Differentiation of pluripotent embryonic stem cells into the neuronal lineage in vitro gives rise to mature inhibitory and excitatory neurons[END_REF] and extend them to our particularly defined culture conditions. Immunostaining and quantitative RT-PCR were used as classical and established methods to link chip-based transcript profiling to other experiments that have been performed with much higher replicate number. Notably, although the shape of PCR and chip profiles was very similar, the induction levels differed sometimes. Such variations may be due to different hybridization efficiencies and background levels. Therefore, very extensive studies, involving RT-PCR controlled by internal standards, will be necessary for a quantitative definition of a final set of markers.

Notably, we did not use differences in absolute numbers of regulations in the present study as basis for any of our conclusions. In addition, all major conclusions are built on groups of coregulated and biologically linked genes as opposed to speculations based on the presence or absence of a single gene with possibly inappropriate oligonucleotide probe. In this context, we are also aware of the issue, that the levels of mRNA do not essentially correlate with the expression, modification and localization of protein, which eventually determines the cellular phenotype. Even if such correlations were given, our approach should not be interpreted as phenotype definition on single cell resolution. This would require different approaches, that are able to measure several markers within one cell [START_REF] Kawaguchi | Single-cell gene profiling defines differential progenitor subclasses in mammalian neurogenesis[END_REF] . The genes grouped within the clusters described here are not necessarily expressed in the same cell (e.g. cluster V: Gfap (astrocytic) and syntaxin (neuronal)) and therefore do not automatically describe a single biological entity. However, with these caveats, we feel that markers that plausibly describe relevant biological processes, and that can be sensitive indicators of disturbances of the default development, can still be selected with confidence on the basis of our study.

In the area of developmental toxicology and especially in DNT, cause-effects relationships are still mostly unknown, and human epidemiological data are only available for a handful of industrial chemicals [START_REF] Grandjean | Developmental neurotoxicity of industrial chemicals[END_REF] . Rodent data based on the OECD test guideline 426 are at least available for around 200 substances [START_REF] Makris | A retrospective performance assessment of the developmental neurotoxicity study in support of OECD test guideline 426[END_REF] . With this lack of human-relevant information and the better animal data base, it appears reasonable to us to perform proof-of-principle experiments for the usefulness of a new approach in rodent cells first, and to validate human cells against these in case of a positive outcome. Moreover, mESC have been shown to model even complex aspects of rodent brain development faithfully, as far as neuronal specification is concerned 2 .

At present, DNT studies are based on e.g. behavioral, cognitive or neuropathological endpoints, and the next step towards mechanistic information would be an understanding of changes on the level of cells and gene expression. The overlapping waves defined here would provide a conceptual framework for this. Such waves (i.e. spatially shifting activation) of gene expression are known from many pioneering studies of mammalian in vivo CNS development [START_REF] Rao | Developmental Neurobiology[END_REF] and are for instance well-characterized in high density and resolution in the hippocampus [START_REF] Mody | Genomewide gene expression profiles of the developing mouse hippocampus[END_REF] . Waves have also been defined in vitro in mESC [START_REF] Aiba | Defining a developmental path to neural fate by global expression profiling of mouse embryonic stem cells and adult neural stem/progenitor cells[END_REF][START_REF] Schulz | The FunGenES database: a genomics resource for mouse embryonic stem cell differentiation[END_REF] or differentiating embryonic carcinoma cells [START_REF] Wei | Global gene expression patterns during neural differentiation of P19 embryonic carcinoma cells[END_REF][START_REF] Przyborski | Transcriptional profiling of neuronal differentiation by human embryonal carcinoma stem cells in vitro[END_REF] . Here, we extended this concept, by relating regulation clusters to underlying biological processes important for toxicity testing. This translation from developmental biology to the toxicological perspective defines the windows of sensitivity relevant for test protocols In the field of cardiac development, the mESC based embryonic stem cell test (EST) has been thoroughly validated and frequently applied [START_REF] Marx-Stoelting | A review of the implementation of the embryonic stem cell test (EST). The report and recommendations of an ECVAM/ReProTect Workshop[END_REF] . Exposure of cells during the entire test period can result in data that are hard to interpret and that are confounded by relatively unspecific toxicity. Therefore separation of exposure into the proliferation and differentiation phase has been suggested [START_REF] Van Dartel | Nervous system development 51 3e-14 Neurod4, Nes (nestin), Cdh2 (N-cadherin), Fgf5, Sema5b, Efnb2 Regulation of nervous system development 17 7e-09 Nefm (neurofilament M), Chrna3 (cholinergic R.), Ntrk3, Isl1, Foxg1 Regulation of neurogenesis 16 9e-09 Hoxa2, Smo, Dll1 (delta-like 1), Hes3, Metrn, Ntrk3 (= Trkc) Neuron projection morphogenesis 13 2e-05 Epha7, Mtap1b, Myh10 (myosin heavy chain), Egr2, Epha7, Isl1 Central nervous system development 21 1e-06 Mtap1b (microtubule-associated protein), Bmi, Foxg1, Isl1, Fgfr3 Neuron projection regeneration 5 2e-06 Mtap1b, Bcl2, Smo, Chst3 (carbohydrate sulfotransferase) Parasympathetic nervous system development 4 4e-06 Hoxb2, Egr2, Smo (smoothened), Hes3 (hairy and enhancer of split) Neuron development 20 5e-06 Mtab1b, Foxg1, Epha7 (Eph receptor A7), Isl1, Ulk2, Bmpr1b Cranial nerve development 5 1e-05 Gli3, Hoxb2, Egr2 (early growth response), Smo, Hes3 Dorsal/ventral pattern formation 9 7e-07 SP8, Foxg1, Bmpr1a, Bmpr1b (bone morphogenic protein R.), Hoxa2 Tissue development 28 1e-06 Homer1, Prox1 (prospero-related homeobox 1), Fzd2, Sdc1 (syndecan) MAPKKK cascade 12 2e-06 Mapk8, Fgf13, Jak2, Nrg1, Fgfr3, Tgfbr[END_REF] . We want to expand this principle here by suggesting four relevant test periods. DoD1-7: testing of lineage commitment, efficiency of NPC formation, axon formation and guidance, and of epigenetic changes associated with the transition from pluripotent cells to more committed NPC. DoD8-15: major phase of neuronal patterning and vesicle development. DoD15-20: a more unexpected, but highly interesting and relevant phase, when most proliferation has ceased, and maturation becomes evident by expression of matrix components, important transporters and disease-associated genes. Our data on lead exposure during this phase show that it will be of high importance for future testing. DoD20+ has not been explored here. It requires, however, further investigation to determine whether this period can be used as stable reference for neurotoxicity vs DNT, or whether new processes such as synaptogenesis, gliogenesis, or myelination take a dominant role here.

Examples for markers of NPC and neuronal differentiation, and the issues associated with such a selection have already been compiled [START_REF] Kuegler | Markers of murine embryonic and neural stem cells, neurons and astrocytes: reference points for developmental neurotoxicity testing[END_REF] . The major task for the future will be the validation of a larger set of such markers, first with known specific and mechanisticallydefined disruptors of developmental pathways, then with known DNT compounds, in order to select the smallest group of final markers useful for a comprehensive description of toxicities triggered by the test compounds.

Materials and Methods

Materials

Cell culture and differentiation

The murine embryonic stem cell (mESC) line CGR8 [START_REF] Suter | A Sox1 to Pax6 switch drives neuroectoderm to radial glia progression during differentiation of mouse embryonic stem cells[END_REF] , kindly provided by K.-H. Krause (Geneva), was cultured in complete Glasgow's modified Eagles medium (GMEM), supplemented with 10% heat inactivated fetal bovine serum (FBS; PAA, Pasching, Austria), 2 mM Glutamax, 100 µM non-essential amino acids, 50 µM β-mercaptoethanol, 2 mM sodium pyruvate and 1000 U/ml leukemia inhibitory factor (Chemicon). Cells were kept at 37°C in 5% CO 2 on tissue culture plates coated with 0.1% gelatin, and were routinely passaged every 48 h.

The mESC were differentiated towards the neural lineage according to the protocol developed by Ying and colleagues [START_REF] Ying | Defined conditions for neural commitment and differentiation[END_REF] . At critical steps, we used the following parameters: cells were plated in the priming phase at 1.2 x 10 5 cells/cm 2 in complete GMEM on 0.1% gelatin coated Nunclon culture dishes (Nunc, Langenselbold, Germany). Next day, for neural induction, cells were plated on gelatin-coated Nunclon dishes at 10 4 cells/cm 2 in N2/B27 medium (composition as described in [START_REF] Ying | Defined conditions for neural commitment and differentiation[END_REF] , for a detailed description of B27 see http://www.paa.com/cell_culture_products/reagents/growthsupplements/neuromix.html). On day 7 of differentiation (DoD7) for neuronal generation and maturation, cells were replated at 10 4 cells/cm 2 on poly-L-ornithin (10 µg/ml) and laminin (10 µg/ml) coated Nunclon dishes in N2/B27 medium. Cells were fed every other day with complete medium change with N2/B27 medium.

Immunostaining and FACS analysis

For immunocytochemical analysis, cells were fixed with methanol (-20°C) or 4% paraformaldehyde (PFA) in phosphate-buffered saline (PBS) and permeabilized with 0.1% Triton X-100 in PBS. After blocking with 10% FBS, cells were incubated with primary antibodies (Tuj1 1:1000, NeuN 1:200, GAD65 1:200, SV2 1:200, PSD95 1:500, Nestin 1:500, Nestin-647 1:40, GFAP 1:800) over night. After incubation with appropriate secondary antibodies nuclei were counterstained with Hoechst H-33342 dye. Images were taken on the original cell culture dishes using an IX81 inverted microscope (Olympus, Hamburg, Germany) equipped with a 40x, NA 0.6 long range lens and processed using CellP imaging software (Olympus). For confocal microscopy cells were grown on 4-well chamber slides (Nunc), fixed with 4% PFA/2% sucrose in PBS and permeabilized with 0.6% Triton X-100 in PBS. After blocking with 5% BSA/0.1% Triton X100 in PBS cells were incubated with Tuj1 antibody in blocking buffer for one hour at room temperature. After incubation with appropriate secondary antibodies nuclei were counterstained with DAPI. Confocal images were taken using a Zeiss LSM 510Meta confocal microscope equipped with a Plan Apochromat 63x, NA 1.4 oil DIC lens. Images were analyzed and processed using ImageJ.

For flow cytometry, cells were dissociated on DoD7 with accutase, fixed and permeabilized in Cytofix Buffer followed by Perm Buffer I (both BD Bioscience, Franklin Lakes, USA), and stained with anti-nestin antibody conjugated to Alexa-647, or isotype control. Cells were analyzed with an Accuri C6 flow cytometer (Accuri Cytometers, Ann Arbor, USA) and data processed with CFlow Plus (Accuri Cytometers).

Quantitative PCR and quality control of differentiation

Total RNA was isolated at indicated time points for marker gene expression analyses using Trizol, the RNA was retro-transcribed with SuperScript II reverse transcriptase, and the resultant cDNAs were amplified in a Biorad Light Cycler (Biorad, München, Germany) with primers specific for the genes of interest and designed for a common melting temperature of 60°C. Real-time quantification for each gene was performed using SybrGreen and expressed relative to the amount of gapdh mRNA using the 2^(-Delta Delta C(T)) method [START_REF] Livak | Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2-[Delta][Delta]CT Method[END_REF] . For each run, the consistency of conditions and constancy of gapdh amounts in the samples was controlled by assessment of its absolute cycle number (= 18 ± 0.5).

Gene expression analysis

Cells were used for RNA preparation as undifferentiated mESC before the priming phase (day 0), on DoD7 (before replating), on DoD15 and on DoD20. RNA was extracted from Trizol preparations and purified using RNeasy Mini prep columns (Qiagen). The total RNA harvested was quantified using a Nanodrop device (Thermo Scientific, USA) and its integrity was assessed using Agilent Bioanalyser (Agilent, USA). Illumina TotalPrep RNA Amplification Kit (Ambion, USA) and 500 ng total RNA of each sample was used according to the manufacturer's protocol to produce biotin-labelled cRNAs.

Data analysis

Original and processed data have been deposited for public access in the EBI Arrayexpress database (Accession Number to be supplied). For initial processing, data were uploaded to Beadstudio (Illumina) for background subtraction. Further processing (baseline transformation and normalization to 75 percentile) and analysis was performed with Genespring 9.0 (Agilent, Santa Clara, CA), and all normalized expression kinetics data sets were used as input for an unsupervised non-hierarchical clustering with relation to the average of expression of all genes on the chip, using the K-means algorithm. The eight major clusters were selected for further analysis. Within these, significant gene expression profiles were selected, based on a minimum regulation of 2.0-fold on any of the time points and on twoway ANOVA taking into account the regulation range and the variation between different arrays.

Patch-clamp recording

Statistics and data mining

The numbers of replicates of each experiment are indicated in figure legends. Data were presented, and statistical differences were tested by ANOVA with post-hoc tests as appropriate, using GraphPad Prism 4.0 (Graphpad Software, La Jolla, USA). Assignment of significantly overrepresented gene ontology (GO) categories to different clusters, and calculation of probabilities of a false-positive assignment was performed by G-profiler (http://biit.cs.ut.ee/gprofiler/ [START_REF] Reimand | Profiler--a webbased toolset for functional profiling of gene lists from large-scale experiments[END_REF] ). For coverage of biological domains without appropriate and well-controlled GO category, relevant genes were assembled from the literature and crosschecked by 2-3 independent specialists. The number of genes within these groups identified in this study was indicated in relation to the overall number of possible hits or in relation to their distribution over different clusters. The genes defined in this study as embryonic stem cell markers or neural stem cell (NPC) markers were derived from recent literature [START_REF] Kuegler | Markers of murine embryonic and neural stem cells, neurons and astrocytes: reference points for developmental neurotoxicity testing[END_REF] . Neuronal (N) differentiation markers (n = 574) were defined as all members of gene ontology (GO) GO:0048699 (generation of neurons) corrected for those genes used as NPC markers. The graphical representation of identified genes (or groups) within their biological context is based on the major gene function as indicated on the NCBI-gene website and the literature.

Importantly, members of each identified group were scored according to their suitability as markers for a PCR-based quality control of the differentiation pattern in toxicity experiments.

Several selection rounds were run to identify the final set of markers displayed as example genes in the tables and some of the figures.

Toxicity experiments

Cells were exposed to chemicals during different phases of differentiation to test the suitability of the model system for neurotoxicity testing, and for testing of developmental neurotoxicity during defined time windows. Retinoic acid (1 µM), "3i" (a mixture of 0.8 µM PD184352, 2 µM SU5402, 3 µM CHIR99021) [START_REF] Ying | The ground state of embryonic stem cell self-renewal[END_REF] or cyclopamine (1 µM) were added to cultures from DoD1-DoD7 or from DoD8-DoD15. Then the experiment was ended, or incubation continued in the absence of chemicals for additional 6 days. On the final day, RNA was prepared by the Trizol method for PCR analysis. For morphological observations, the monolayer regions within the culture wells were imaged. Genes were preselected before the analysis as endpoints for initial proof-of-concept experiments, and results from all genes chosen are presented. A. Gene expression kinetics were determined for all genes represented on the chip. An unbiased clustering analysis of the kinetic profiles of all regulated genes was performed. For each cluster (named Ia, Ib, IIa, IIb, IIIa, IIIb, IV, V), the means of the absolute expression level of all genes in the respective cluster, for each analysis time point is displayed. and plotted on a logarithmic scale; n: number of genes in the cluster B. Number of genes expressed in mESC, NPC and developing neurons (N) were analyzed by extensive literature search (mESC, NPC) or GO-analysis (N). The relative distribution of these genes across the different clusters was calculated (in %) and displayed (e.g. 65% of all ESC markers were found in cluster Ia, 35% of all N markers in cluster III).

  Unless otherwise mentioned, cell culture media and reagents were from Invitrogen (Carlsbad, USA) and accessory reagents from Sigma. Antibodies: anti-Tuj1 (cat. # MMS-435P; Covance), anti-NeuN (cat. # MAB377; Chemicon), anti-GAD65 (GAD-6; DSHB), anti-SV2 (SV2; DSHB), anti-PSD95 (cat. # 51-6900; Zymed), anti-Nestin (cat. # MAB353; Chemicon), anti-GFAP (clone: G-A-5; Sigma) anti-Nestin-647 (clone: 25/NESTIN; BD Biosciences).

  mRNA Primer: Pou5f1-fw: CTCTTTGGAAAGGTGTTCAGCCAGAC, Pou5f1-re: CGGTTCTCAATGCTAGTTC-GCTTTCTC; Nestin-fw: CTGGAAGGTGGGCAGCAACT, Nestin-re: ATTAGGCAAGGGGGAAGAGAAGGATG; Synaptophysin-fw: GGGTCTTTGCCATCTTCGCCTTTG, Synaptophysin-re: CGAGGAGGAGTAGTCACCAACTAGGA; Gfap-fw: GCCCGGCTCGAGGTCGAG, Gfap-re: GTCTATACGCAGCCAGGTTGTTCTCT; Shh-fw: CAGCGGCAG-ATATGAAGGGAAGATCA, Shh-re: GTCTTTGCACCTCTGAGTCATCAGC; Hes5-fw: CCCAAGGAGAAAAACC-GACTGCG, Hes5-re: CAGCAAAGCCTTCGCCGC; Tubb3: GACAACTTTATCTTTGGTCAGAGTGGTGCTG, Tubb3re: GATGCGGTCGGGGTACTCC; Nkx2.1-fw: TACCACATGACGGCGGCG, Nkx2.1-re: ATGAAGCGGGA-GATGGCGG; Dlx1-fw: TCACACAGACGCAGGTCAAGATATGG, Dlx1-re: AGATGAGGAGTTCGGATTCCAGCC; -re: AGGCGGCTCATTCTCTCTTCATTGT; Isl1-fw: ACCTTGCGGACCTGCTATGC, Isl1-re: CCTGGATATTAGTTTTGTCGTTGGGTTGC, Tubb3-fw: GACAACTTTATCTTTGGTCAGAGTGGTGCTG; Tubb3-re: GATGCGGTCGGGGTACTCC, Mapt-fw: ACACCCCGAACCAGGAGGA; Mapt-re: GCGTTGGAC GTGCCCTTCT ; App-fw: TCAGTGAGCCCAGAATCAGCTACG, App-re: GTCAGCCCAGAACCTGGTCG, Pink1-fw:

  For hybridization onto Sentrix Mouse Ref.8 V2 mRNA microarray beadchips (Illumina), 750 ng labelled cRNA were incubated for 16 h at 58°C. After hybridization, chips were washed, blocked, and streptavadin-Cy3 stained. Fluorescence emission by Cy3 was quantitatively detected using BeadArray Reader Scan. Statistical analysis data is based on duplicate samples. Each of the samples contained pooled RNA from two differentiations to further increase robustness of results. Technical variation of the chip was minimal as tested by rerun of the same sample on two different arrays and by comparison of results from two beadchips within one array.

  Electrodes with a resistance of 2-5 MΩ were pulled of borosilicate glass (Clark, G150F, Warner Instruments, Hamden, CT, USA) on a Sutter Instruments (Novato, CA, USA) P-97 horizontal micropipette puller. All experiments were carried out using a custom built recording chamber (800 µl volume) made of Teflon within a temperature-controlled microscope stage (37°C). Whole cell voltage and current clamp recordings were obtained from cells at day of differentiation (DoD) 20-24. Cells were grown on coated glass cover slips (10 mm) from DoD7 on. Whole-cell currents were recorded using an L/M-EPC-7 amplifier (List Medical Electronic, Darmstadt, Germany), digitised at sampling frequencies between 10 kHz to 50 kHz using a DigiData 1320A AD/DA converter (Axon Instruments Inc.). The patch pipettes for spontaneous and evoked action potential measurements as well as for the neurotransmitter responses were filled with (in mM) 90 K + -gluconate, 40 KCl, 1 MgCl 2 , 10 NaCl, 10 EGTA, 4 Mg-ATP, 10 HEPES/KOH (pH 7.4 at 37°C), whereas the bath solution contained (in mM): 155 NaCl, 1 CaCl 2 , 3 KCl, 10 D-(+)-glucose, 10 HEPES/NaOH (pH 7.4 at 37°C). The protocol for recording of Na + and K + channels was as follows: cells were hyperpolarized to -90 mV, and subsequently stepped to a defined voltage as indicated and returned to -70 mV, before the next cycle a different voltage step was run. Each cycle took 120 ms. For the neurotransmitter response measurements, the different substances were directly added as concentrated stock solutions to the recording chamber in amounts of 1-10 µL. Antagonists were added at least one min before the agonists. Recordings were initiated within 100 ms after addition of agonists. For the measurement of barium currents through calcium channels the pipette filling solution contained (in mM) 110 CsF, 10 NaCl, 20 TEA-Cl, 10 EGTA, 4 Na 2 -ATP, 10 HEPES/CsOH (pH 7.4 at 37°C), whereas the bath solution contained (in mM) 130 NaCl, 10 BaCl 2 , 10 D-(+)-glucose, 5 tetraethylammonium chloride, 10 4-aminopyridine, 0.5 tetrodotoxin, 10 HEPES/NaOH (pH 7.4 at 37°C).All current signals were normalized against the individual cell capacitances (as a surrogate measure for cell size) and are expressed in current densities (current divided by cell capacitance). Liquid junction potentials (LJP) were measured and corrected, using the method described by Erwin Neher (1992) except for barium current measurements. Stimulation, acquisition and data analysis were carried out using pCLAMP 10.2 (Axon Instruments Inc.) and ORIGIN 8.0 (OriginLab Corp., MA, USA). Fast and slow capacitive transients were cancelled online by means of analogue circuitry. Residual capacitive and leakage currents were removed online by the P/4 method. Series Resistance Compensation was set to at least 50%. For analysis, traces were filtered offline at 5 kHz. Cells for measurements were chosen with respect to their morphological phenotype (small round highly elevated (phase-bright) cell bodies with projections of at least five times cell body diameter, growing in network-like clusters containing at least 20-30 similar cells). The patch pipette was approached to these cells perpendicular to the plane formed by the cell membrane in the patch region.

Figure 1 .

 1 Figure 1. Protein and mRNA-based markers of robust neuronal differentiation

Figure 2 .

 2 Figure 2. Electrophysiological evidence for successful neuronal development.

Figure 3 .

 3 Figure 3. Detection of non-cytotoxic developmental disturbances by

Figure 4 .

 4 Figure 4. Cluster analysis of mRNA time course profiles, and their association

Figure 5 .

 5 Figure 5. Indication of a progressive change in chromatin organization and

Figure 6 .

 6 Figure 6. Correlation of neural precursor formation with a transiently

Figure 7 .

 7 Figure 7. Analysis of glia-associated genes

Figure 8 .

 8 Figure 8. Functional assignment of neuronal genes up-regulated in different

Table 1 :

 1 GO categories significantly overrepresented in cluster IIa* All categories identified by gProfiler bioinformatics analysis, with their p-values indicated after correction by removal of "nervous system development" genes from non-neuronal GOs.

	Biological process (GO)*	genes	p-
		in IIa	value

examples of upregulated genes listed in the GO

Acknowledgement

The Work was supported in part by the Doerenkamp-Zbinden Foundation, the DFG, the EU FP7 project ESNATS (ML, SK), an IRTG1331 fellowship (BZ) and a fellowship from the KoRS-CB (PBK). We are grateful to Giovanni Galizia and Sabine Kreissl for help with the electrophysiological recordings and indebted to Bettina Schimmelpfennig for invaluable experimental support. The monoclonal antibodies Gad-6 developed by D.I. Gottlieb and SV2 developed by K.M. Buckley were obtained from the Developmental Studies Hybridoma Bank developed under the auspices of the NICHD and maintained by The University of Iowa, Department of Biology, Iowa City, IA 52242. We thank K.H. Krause for the CGR8 mESCline and J. Vilo and S. Ilmjärv for help with bioinformatics analysis.

Conflict of interest

The authors declare no conflict of interest


Eui-Jong Kim

Jean-Jacques Roux

Gilles Rusaouen

Frédéric Kuznik
email: frederic.kuznik@insa-lyon.fr

Numerical modelling of geothermal vertical heat exchangers for the short time analysis using the state model size reduction technique

Keywords: borehole heat exchanger, numerical simulation, state model size reduction, short time analysis, FEM

According to the introduction of a dynamic operating mode in ground-coupled heat pump systems, a short time analysis within and around borehole heat exchangers is required in the modern geothermal system simulation. A numerical modelling could be a proper answer for this challenge. However, the numerical model is time consuming and necessitates a large memory particularly in such large systems. Therefore, the state model size reduction technique has been applied in this paper with various numerical techniques particularly in the finite elements method. As a result, the reduced model developed is: a) relevant with a validation using a traditional analytical model (using 100% modes) and b) efficient in calculation time, only using 6% modes and consequently reducing time consumption up to 95%.

Introduction

Since many years, the installation of ground-coupled heat pump (GCHP) in buildings has increased because of their efficiency for heating and cooling. Contrary to this aspect, the relatively higher initial cost of geothermal vertical heat exchangers, or called borehole heat exchangers (BHE) would have been an obstacle to the spread of the systems particularly in the residential buildings.

On the other hand, a combination of the GCHP with solar thermal collectors has been studied in recent years, for example the direct thermal recharge of the soil through the same U tubes using the solar collector [START_REF] Trillat-Berdal | Experimental study of a ground-coupled heat pump combined with thermal solar collectors[END_REF]. Chiasson [START_REF] Chiasson | Assessment of the viability of hybrid geothermal heat pump systems with solar thermal collectors[END_REF] showed in his simulation analysis that the combination system makes possible to reduce BHE length. A dynamic operating mode being introduced in such system, a short time analysis within and around the borehole is required for understanding heat transfer having characteristic time of one hour or even few minutes.

In order to give an appropriate response to this challenge, several models have already been developed, based on analytical or numerical solutions. Lamarche et al. [START_REF] Lamarche | New solutions for the short-time analysis of geothermal vertical boreholes[END_REF] proposed a model which is capable to analytically take into account the effect of inertia of the borehole for the short time analysis. In spite of their exactness under the given conditions, the developed analytical models cannot be easily used for any conditions. Moreover, many coefficients used in the analytical model can only be obtained a posteriori using a numerical model or experimental data. Despite everything, numerical analyses for this purpose have been continuously developed, for example the work of Yavuzturk et al. [START_REF] Yavuzturk | A Transient Two-Dimensional Finite Volume Model for the Simulation of Vertical U-Tube Ground Heat Exchangers[END_REF] which was considered as the one of the first numerical approach modelling geothermal vertical boreholes, i.e. the pipes and the grout with their inertia. More recently, Cui [START_REF] Cui | Numerical analysis and experimental validation of heat transfer in ground heat exchangers in alternative operation modes[END_REF] compared his FEM (Finite Element Method) to the traditional analytical model, in the case of the infinite line source [START_REF] Ingersoll | Heat Conduction with Engineering, Geological, and Other Applications[END_REF].

The results prove the importance of the short time analysis particularly in an alternative operating schedule. Although the FEM numerical method is accurate, the relatively large calculation time and the large memory required are its shortcomings for an application in such large system especially in the case of a long period simulation (typically a year for building application after many-year-operation of the system).

Therefore, in this paper, we propose a numerical model for short time transient response, using a state model size reduction technique for decreasing the calculation time. The part 2 of the paper details the model size reduction technique. The numerical model is presented in the third part of the article and the fourth part deals with the modelling strategies. The part 5 presents the results of the numerical simulations showing the method efficiency.

State model size reduction technique

For thermal analysis, model size reduction techniques have been used for different purposes over the last 30 years since its invention [START_REF] Marshall | An approximate method for reducing the order of a linear system[END_REF] in the field of automatic controls. Among them, several transient analyses of large systems found in heavy building components have been studied. For example, Mokhtari [START_REF] Mokhtari | Etude du comportement thermique des bâtiments en contact avec le sol par l'analyse modale[END_REF] studied the thermal interaction between the building floor and the surrounding soil; Gao et al. [START_REF] Gao | Dynamical building simulation: A low order model for thermal bridges losses[END_REF] proposed a model for the thermal bridge in the wall corner with a 3D approach. Then, we will use the model size reduction technique whose implementation process will be described in the following parts.

State space model

Above all, some assumptions are supposed for our mathematical modelling as follows:

-The soil property ( , , c) is considered constant; consequently, constant -The heat transfer is considered 2-dimensional in the Cartesian -The perfect thermal contact between the materials is supposed Therefore, the heat equation can be expressed as:

2 2 0 0 0 2 2 ( , , ) ( ) T T T x y t t x y α ∂ ∂ ∂ = + ∂ ∂ ∂ (1)
For numerical analysis, the continuous problem must be converted into an algebraic problem.

Then, the spatial discretization of Equation 1 using finite volume method (FVM) or finite element method (FEM) leads to the following linear equation system:

) ( ' ) ( ) ( 1 1 t U B C t AT C t T - - + = (2)
Where T(t), the approximated state vector (dimension: n, 1), is the assembly of the time dependant nodes temperature of the system, that is, T(t) contains all the nodes temperatures which locations are defined by the spatial discretization :

T(t)={T 1 (t), T 2 (t), T 3 (t), … T n (t)}',
where 1,2,3…n indicate the node number order in the mesh generation. The matrix C, the square capacitance matrix (dimension: n, n), is related to the thermal inertia of the finite volumes, and the matrix A, the square conductance matrix, contains thermal relations among the nodes represented by the state vector T(t). U(t), the solicitation vector (dimension: p, 1), represents the boundary conditions of the system, therefore; the matrix B , a rectangular matrix (dimension: n, p), represents physical relation between the boundary conditions and the system.

In dynamic simulations, only the evolution of some particular values is necessary: for example, temperatures or heat flux in specific boundaries. For these values, another state model could be defined as:

) ( ) ( ) ( t DU t JT t Y + = (3)
The output vector Y(t) (dimension: q, 1) is the assembly of the concerning outputs, J, the observation matrix (dimension: q, n) and D, the direct transmission matrix (dimension: q, p) are defined under the physical correlation with the state vector and the solicitation vector, respectively.

The state model mentioned above is a complete expression for the dynamic simulation.

However, the numerical approach for the fine analysis of such large system is always calculation time consuming. Therefore, the state model size reduction technique is naturally expected, and then the reduction processes is be described in the following part.

Reduction technique

The objective of the state model size reduction technique is to reduce the order of the abovedescribed state space model (m << n) (see equation 2-3) and then the calculation time. There are several approaches already developed for various purposes, so we will only present the method used in this paper, demonstrating a common philosophy concerning the category (modal basis truncation) to which our method belongs. The different categorizations were proposed by Petit [START_REF] Petit | Réduction de modèles de connaissance et identification de modèles d'ordre réduit[END_REF].

Above all, a basis change into the modal basis has to be accomplished ahead of the reduction processes because the new model reproduced in the modal basis could be truncated. The modal basis model could be obtained by the diagonalization of the matrix C -1 A (2).

Consequently, the modal basis model and its correlations from the equations ( 4) and ( 5) can be written as:

+ Ω = + = ) ( ) ( ) ( ) ( ) ( ) ( t DU t X t Y t BU t WX t X (4) JP B C P B AP C P W t PX t T = Ω = = = - - - - , ' , ), ( ) (
eigenvalue i of the initial matrix C -1 A: W=diag( i ), and the associated eigenvectors V i compose the matrix

P: P=[ V 1 V 2 …..V i ….V n ].
Then, various reduction methods for Equation 4 can be used. Marshall [START_REF] Marshall | An approximate method for reducing the order of a linear system[END_REF] firstly proposed a method for choosing the dominant modes, i.e. the critical modes for deciding the dynamic evolution of the outputs (Y(t)). He introduced a temporal criterion in order to select the modes that are the most meaningful for the dynamic analysis: i = i -1 > t/4 with t the time step used in the numerical simulation. The rest modes could be used only for the constitution of the steady state condition. This method is simple but inefficient for the common purpose of truncating most of the non-dominant modes. Another method often used in the building thermal analysis is the linear aggregation method. This method uses the notion of energy contribution whose matrix terms are defined by the relation [START_REF] Ingersoll | Heat Conduction with Engineering, Geological, and Other Applications[END_REF]. The criterion for selecting the dominant modes, i.e. the most energetic modes, is:

= = ≅ n i ii m i ii E E 1 1
with m the number of the dominant modes. The details about the construction of E ij are given in [START_REF] Michaïlesco | Approximation des systèmes complexes par des modèles de dimension réduite[END_REF], where r has 0 for Dirac impulse or 1 for Heaviside step function. In the truncation category, the method of Moore was developed by a slightly different approach that is detailed in [START_REF] Moore | Principal component analysis in linear system: controllability, observability, and model reduction[END_REF].

) ( ) ( 1 1 j i r j i q l p k jk ik lj li ij B B E ξ ξ ξ ξ + Ω Ω - = ′ = = (6)
In the present study, the reduction uses the linear aggregation method but the final equations [START_REF] Marshall | An approximate method for reducing the order of a linear system[END_REF] would have the same structure for all the techniques briefly presented previously.

First, the unselected modes are removed from the initial equation ( 4) and, consequently, the state order can be reduced from n to m with m<< n. Therefore, Y(t) could be obtained with approximated values (t) but with faster calculation. The error due to this reduction processes could vary according to the chosen method and the error definition in the criterion... As we already mentioned, the non-dominant modes unselected just contribute to the accuracy of the calculated steady state (see the last term in each equation of ( 7)).

Ω - + Ω = - = + = - - - - - - - - ) ( ) ( ) ( ) ( ~) ( ), ( ) ( ) ( 1 1 t U G W D t X t Y t U G W X t U G t X W t X m n m n m n m m m n m n m n m m m m (7)

Numerical model

To implement the state model size reduction technique, a specific computer code has been written in spite of the availability of commercial simulation packages because the reduction procedure contains many matrix transformations. In the following parts, the numerical methods employed for the new code will be described.

Unstructured mesh

As mentioned in the introduction, a fine mesh should be used in this model for the exact short time response of the geothermal vertical boreholes. Many various meshes have been proposed like the structured mesh [START_REF] Yavuzturk | A Transient Two-Dimensional Finite Volume Model for the Simulation of Vertical U-Tube Ground Heat Exchangers[END_REF], [START_REF] Rottmayer | Simulation of a Single Vertical U-Tube Ground Heat Exchanger in an Infinite Medium[END_REF] or the unstructured mesh [START_REF] Cui | Numerical analysis and experimental validation of heat transfer in ground heat exchangers in alternative operation modes[END_REF], [START_REF] Li | Development of a numerical model for the simulation of vertical U-tube ground heat exchangers[END_REF]. However, the structured mesh models could not describe the circular configuration found within the boreholes. For example, Rottmayer [START_REF] Rottmayer | Simulation of a Single Vertical U-Tube Ground Heat Exchanger in an Infinite Medium[END_REF] proposed a geometric factor for reducing the errors caused by adopting the equivalent non circular shape. Moreover, the structured models developed in the polar grid system cannot be flexible for a multiple borehole simulation.

For these reasons, an unstructured mesh based on the Delaunay triangulation for the horizontal plan was adopted in the paper (see Figure 1). In this work, the mesh is generated by Gmsh [15]. The mesh generator outputs mesh information (i.e. a mesh file) so that the grid information could be integrated in our code with an additional identification procedure.

Finite elements method

On the one hand, the unstructured mesh is a good solution because of the flexible mesh generation and consequently the possibility to obtain a lower mesh number in a complex configuration system. On the other hand, the non-orthogonal relationship caused by this unstructured mesh must be corrected, particularly when using FVM, which has been a common method in the conduction heat transfer analysis. However, the correction method, using, for example, the least-square approach [START_REF] Demirdži | Numerical method for coupled fluid flow, heat transfer and stress analysis using unstructured moving meshes with cells of arbitrary topology[END_REF], could make the conductance matrix A (Equation 2) asymmetric. In addition, an approximation with FVM like average thermophysical properties between two different materials must be used in the correction formula although a special correction method would have been used for the symmetry of A. As a result, the FEM, a method free from these approximations by its nature, is applied in the model.

For FEM, the weak formulation of the governing equation (1) could be achieved by using some mathematical modification using a weighting function T*. The obtained low-order Equation 8 is not handling the non orthogonal nature in the process of the matrices definition.

* * * 0 dV dV T T T T t c T V V ∀ = ∇ ∇ + ∂ ∂ λ ρ (8)
Ahead of transforming Equation 8 into the algebraic equation (Equation 2), an approximation of temperatures inside an element is supposed by an interpolation function (Equation 9-10). In this paper, all the expressions are described for the 2D triangulation mesh.

{ }

1 1 2 3 2 3 ( ) ( , , ) ( , ) ( , ) ( , ) ( ) ( , ) ( ) ( )
∆ + + = 2 ) , ( y c x b a y x N i i i i (10)
Where, a i , b i , c i will be defined by the geometric information in each element, and the area of the element considered.

For the temperature gradient, another function could be employed from N i formulation.

∆ = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = 3 2 1 3 2 1 3 2 1 3 2 1 2 1 ] [ c c c b b b y N y N y N x N x N x N M ( 11
)
From these approximations, the weak formulation (Equation 8) for the element could be expressed as Equation 12, adding a flux (W/m) boundary condition on pipe walls if the element belongs to a pipe wall mesh.

{ } { } { } { } 0 ) ] ' [ ] [ ] [ ( * = - -
∆ ∆ = d N N c C T e ] [ ] [ ] [ ρ (13) ∆ ∆ - = d M M A T e ] ][ [ ] [ ] [ λ (14) dl N B l T e - = ] [ ] ' [ (15)
For the triangulation mesh, the matrices (Equation 13-15) can also be expressed as:

∆ = 2 1 1 1 2 1 1 1 2 12 ] [ c C e ρ (16)
+ + + + + + + + + ∆ - = 2 3 2 3 2 3 2 3 1 3 1 3 3 2 3 2 2 2 2 2 1 2 1 2 3 1 3 1 2 1 2 1 2 1 2 1 4 ] [ c b c c b b c c b b c c b b c b c c b b c c b b c c b b c b A e λ (17) = 1 1 2 ] ' [ l B e (18)
All these matrices (Equation 16-18) defined for each element in the considered domain have to be assembled by an element-by-element overlap (Equation 19) for completing the algebraic equation (Equation 2).

Modelling strategies

A complete BHE model could be developed with the methods defined in the previous sections. However, a naïve numerical approach is not desirable in terms of the calculation time and the occupying memory because of the presence of large size matrices. In fact, the required calculation time for an n-order matrix is generally in the order of n 3 and the necessary memory is increased by the same order. Accordingly, some modelling strategies are necessary.

These strategies will be developed in the followings.

Slice model

In this paper, a novel slice model of the ground soil volume is introduced and accomplished by extruding a horizontal cross-section as shown in Fig 2 . This slice model could be used for a complete model of the relatively long borehole simulation (case1): when t H 2 /100 [START_REF] Eskilson | Thermal analysis of heat extraction boreholes[END_REF]. In addition, several uses of this pre-defining slice model are suitable for describing the fluid temperature variation with borehole depth (case 2). For the 3D simulation, accounting for the thermal interaction among the slice models, the model could be employed just by adding the temperatures of the adjacent models to the solicitation vectors U (case 3). This strategy was motivated by the fact that each slice model has normally an identical mesh configuration, so that the definition of the corresponding matrices (Equations 2-5 and 7) for one slice model is sufficient for its multiple use and particularly for the reduced state model. In addition, the calculation time and the necessary memory are dramatically reduced if partial models are employed for the entire model: ( )

3 3 a Nn bNn >>
where N is the number of slice models and n, the number of nodes in a slice. In this article, only the 2-dimensional slice model (case 1 in Fig. 1) is treated as a complete model, and then this model results would be applied directly for the other applications (case 2-3).

Decomposition of the slice in sub-zones

Following the conclusion of the previous section, the model efficiency in calculation time is enhanced with the sub-zone decomposition of the model presented in Figure 3. In this paper, 4 sub-zone models were proposed, which cover packed soil regions respectively (see Figure 3-b). This method makes it possible to use different time steps for each sub-zone, which could generally be helpful for the time reduction. Additionally, this decomposition is indispensable to take advantage of the state model size reduction technique (see the section 5). The decomposition has also the benefit that, in the presence of multiple boreholes, the matrices definition and their modification for a borehole sub-zone could be reused for other boreholes sub-zones (see Figure 3-c).

In order to decompose the model, a coupling method must be applied to interconnect physical frontier values like temperatures and fluxes until each value converges in the given time step and for the domain considered. However, the nodes on the frontier mesh in the decomposed model have some information of the previously neighboring elements which are currently in the adjacent sub-zone model. Therefore, a special attention must be paid to conserve all the elements information, particularly in the nodes on a decomposition frontier when using decomposition in FEM. In brief, a new temperature boundary has been defined between the sub-zone models decomposed (see Figure 4). With the boundary, all elements information around the frontier nodes could be recovered with some matrices modifications. Then, these sub-zone models employ their own time steps, and then the boundary results in each sub-zone time step are interpolated or extrapolated by a linear function in order to well correspond with the adjacent sub-zones time step in the temporal coupling under the acceptable error criterion, and then keep the accuracy of the numerical scheme. Finally, all the spatial and temporal coupling is completed by an iterative calculation.

Boundary conditions

For the boundary conditions, flux is imposed at the exterior walls of pipes, and the upper and lower boundaries of the slice model are adiabatic (case 1 in Figure 2). The fluxes on each tube, each branch of U-tube, could be calculated iteratively by a coupling with a fluid model which will be not treated in this study but replaced by an imposed flux condition. The parameters and boundary conditions used for the simulations are presented in Table 1 and 2 respectively.

Two cases of heat injection are tested: constant flux case for the steady state regime and long time analysis; sinusoidal flux case for the alternative operating analysis (Table 2). For the model size, a 100m radius soil domain case is treated which is suitable for a long time simulation of about 20 years, according to the criterion of Eskilson [START_REF] Eskilson | Thermal analysis of heat extraction boreholes[END_REF]: r max =3( t max ) 0.5 .

Results and discussion

Mesh independency

The mesh independency concerning simulation results is important for numerical solution because it will decide model relevance and lightness in terms of error and calculation time respectively. For this purpose, an error evaluation has been defined, based on the maximum temperature difference between the two compared models, and this evaluation (Equation 20) will be used for other simulation results in this paper.

test ref T T max error Maximum - = (20)
Figure 5 shows the result on the comparison of different mesh models for the mean borehole wall temperature when injecting constant heat. The first three mesh models (n=1918, 1406, 1018) have the maximum error less than 0.1 (°C), correspondent relative error 0.43%,among them, which was comparatively smaller than other models (n=625 and 313). Therefore, the mesh model whose order is n=1018 was selected as the reference simulation mesh model in this article. This mesh configuration is presented in Figure 6 and the correspondent mesh information in Table 3.

Model validation

In order to validate our approach and show the efficiency of the method proposed in this paper, the numerical model is compared to the analytical model of the infinite line source [START_REF] Ingersoll | Heat Conduction with Engineering, Geological, and Other Applications[END_REF].

This test case has been chose because small errors for each time step in a numerical model could normally accumulate for long time analysis in the geothermal system. The figure 7 shows the comparison between the analytical model and the numerical data which are in good agreement.

For detailed information, the relative error [21] is used to estimate the error between the two models:

Numerical value -Analytical value Relative error 100% Numerical value -Initial value = × (21)

The figure 8 shows the evolution of the relative error. Some differences between the two models exist during the first period (Fo= t/r b ²<20) but are normal because the infinite line source model is only valid if the time criterion Fo>20 is satisfied [START_REF] Ingersoll | Heat Conduction with Engineering, Geological, and Other Applications[END_REF].

When using this entire slice model, the reduction technique is not suitable because the spectra of eigenvalues in this case are very large. Then, the slice decomposition into multiple subzone models is mandatory in order to dense the spectra of eigenvalues which normally decide the system evolution characteristic in a dynamic analysis.

Before examining the sub-zone reduced model, the sub-zone models using the coupling method (spatial and temporal) has to be compared to the entire slice model. Fig 9 shows the temperature evolution on the borehole wall for the models examined. The results of the models are in good agreement. The maximum error of the spatial coupling was 0.00038 (°C), and 0.078 (°C) for the spatial and temporal coupling.

Comparison between complete model and reduced model

For the reduced model (RM) simulation, 4 sub-zones are used and each sub-zone has a different simulation time step: 1 hour, 6 hour, 1 day, and 1 week from the borehole to exterior respectively while the complete model (CM) has only one time step (1hr) for all the sub-zones.

The complete model (CM) uses all of the system modes and is compared to the reduced model (RM) using the linear aggregation (LA) technique with 99% energy criterion:

[%]) 100 / (

1 1 × = = = n i ii m i ii E E .
As shown in Figure 10, RM-LA model and CM are in good agreement with the maximum error of 0.31°C. The selected modes ratio of RM-LA model to CM is about 6% (66<<1018)

and, consequently the reduction in computational time is about 90%. This efficiency could be explained by the result presented in figure 11. Each sub-zone model has an unambiguous energy distribution; that is, some primary modes might decide the dynamic evolution in the given time step, which did not appear in the original entire slice model.

Subsequently, a dynamic heat injection is studied for evaluating the validity of the modelling in the case of short time operation. The heat flux is supposed to be sinusoidal as explained in the boundary condition section. The figure 12 shows the pipe wall mean temperature evolution for CM and RM-LA models with different energy criteria in order to observe the energy criterion impact. Above all, the two models (CM and RM-LA 99%) are in good agreement with a maximum error 0.274 (°C). However, other models with 95 and 85% energy criteria have 1.29 and 2.40 (°C) in maximum error respectively. Although the amount of the selected modes decreases according to the energy criterion, the amount of the selected modes is not greatly different among the reduced models. Furthermore, the CPU time for the transient calculation among the different reduced models was almost similar: in this case, 70% calculation time reduction is obtained compared to CM because of convergence problem in such a dynamic input. Therefore, using RM-LA 99% could be a good solution in this decomposed slice model.

In order to further examine the RM precision compared to the CM and generalize the proven results, we use the bode diagrams. For this simulation, the output mean temperature on the pipe wall is observed when the input flux on the same wall is introduced in the borehole region. First of all, the result (see Figure 13) shows that the magnitude difference between CM and RM is almost zero for frequencies lower than 2×10 -2 which corresponds about 5min excitation period. In terms of the phase difference, the two models have nearly the same response in the area lower than the frequency 2×10 -3 . Although this frequency matches about 50min excitation period, this period level would be sufficient for application within geothermal system simulation field even for the short time analysis. Actually, a simulation with a time step 1min with 10 min periodic input fluxes (between the frequencies 2×10 -2 and 2×10 -3 ) showed a good concordance between the CM and RM with a slight error less than the result of Figure 12. Additionally, when using one minute time step, the RM had the advantage that it reduces 95% of the transient calculation time compared to the CM. This efficiency gives a merit to the model developed (RM-LA) particularly in the short time analysis because this efficiency has decreased when the larger time step was applied (70% time reduction in 1 hour time step).

Conclusion

A short time analysis is required in modern BHE simulation by the presence of dynamic operation modes in the order of an hour or even few minutes. For such an analysis it is necessary to take into account the inertia within and just around a borehole for describing the precise transient phenomena. In this context, a numerical model has been developed which could simulate the borehole configuration and would be calculation time efficient. Then, an unstructured mesh generator (Gmsh) has been used and various techniques have been employed: from a slice model via sub-zone structure to the use of reduced models.

Particularly, the finite elements method and the use of different time steps in each sub-zone have been applied to implement the reduced technique. Moreover, a coupling method for FEM has been applied.

The model has been validated using analytical solutions from the literature. The reduction in computational time in the reduced model was up to 95% compared to the complete model.

The main interest of the developed model is that short time dynamic effects can be easily and rapidly simulated, allowing to significantly enhance the understanding of coupled boreholes systems.

The results can also be a complete reference for the short time analysis when coupling a fluid model for the sufficiently long vertical borehole. Another 3D model for the relatively short borehole is under development based on the slice model by multiplying the slices and adding some solicitation vectors to the slice matrices with newly defined boundary conditions like an ambient temperature on the soil surface.

  function T* is an arbitrary function, only the right part of the equation 12 is required to be zero. The capacity matrix C, the conductance matrix A, and the matrix B for the element considered could be expressed by a combination of the approximation functions N and M as:

Figure 1 .Figure 2 . 4 Figure 3 .Figure 5 .Figure 7 .

 124357 Figure 1. Example of the Delaunay triangular mesh by Gmsh

Figure 9 .Figure 10 .Figure 11 .Figure 12 .Figure 13 .

 910111213 Figure 9. Spatial and temporal coupling validation in sub-zone structured model

Table 1 .

 1 Simulation parameters

	Parameters	Unit	value
	Borehole radius	m	0.055
	External U pipe radius	m	0.016
	U pipe distance (centre to centre)	m	0.06
	Pipe number in a borehole		2
	Soil conductivity	W/m°C 3.5
	Grout conductivity	W/m°C 1.3
	Soil diffusivity	m²/s	1.62×10 -6
	Grout diffusivity	m²/s	3.33×10 -7
	Basic time step	hr	1
	Error criterion for iteration	%	10 -4

Table 2 .

 2 Boundary conditions

	Boundaries	Unit	Value
	Initial soil temperature	°C	10
	Far field temperature	°C	10
	Constant flux density in each pipe	W/m² 300
	Sinusoidal flux density in each pipe W/m² 200+100sin(2 t/24hr)

Table 3 .

 3 Mesh information used for the slice model

	Mesh information	Unit	Value
	Total nodes number		1018
	Nodes number in 4 sub-zone models
			474/140/203/201
	(from the borehole zone to the exterior zone)
	Average mesh size	m²	1.5798×10
	Maximum mesh size	m²	3.1906×10 2
	Minimum mesh size	m²	4.1557×10 -6

1 1 1 (5)The modal state X(t) indicates an assembly of modes of the system which related the dynamic evolution of the system, and each element of the diagonal matrix W corresponds to an

Acknowledgements

The authors want to thank EDF (Electricity of France) R & D for its financial, scientific, and technical support in the framework of the French support program CIFRE.

Sub-zone Model