in addition to classical endpoints , toxicogenomics methods have 
 been used to characterize the assays and to classify toxicants . for 
 regulatory purposes , the descriptive data from , e.g. , transcriptomics 
 or metabolomics approaches need to be converted to quantifiable measures 
 allowing one to compare and predict the hazard of chemicals and drugs .
first attempts to determine benchmark 
 concentrations on the basis of transcriptomics data have already been 
 undertaken in vivo.(6,7 ) their application 
 to in vitro test systems is expected to yield important 
 information about low - dose toxicant effects .
the developing 
 central nervous system is one of the most frequent 
 targets of systemic toxicity .
moreover , testing of 
 nervous system development and possible long - term effects is particularly 
 challenging .
animal testing for example according to oecd guidelines 
 414 ( 2-generation reproduction ) or 426 ( developmental neurotoxicity ) 
 is time - consuming and labor intensive .
the testing capacities are 
 by far not sufficient to test all compounds in vivo that should be studied for reproductive toxicity or developmental 
 neurotoxicity ( dnt ) in the context of reach . moreover
, toxicological 
 risk assessment needs to be adapted to the modern needs of the pharmaceutical 
 as well as the chemical industry . between 
 1999 and 2011
, 19% of the 279 new approved drugs in europe were reported 
 to have postapproval safety - issues , and 5 drugs were withdrawn from 
 the market .
the chemical industry , however , 
 is confronted with the reach initiative asking them to provide more 
 detailed toxicological data regarding new compounds as well as on 
 chemicals already on the market .
consequently the national research 
 council of the usa has published their vision of toxicology for the 
 21st century in 2007 in which they favor the development of mostly human 
 based in vitro test systems , which can be used for 
 high - throughput screening and for high content endpoints such as for 
 transcriptomics approaches .
cultures of differentiating pluripotent 
 stem cells , such as human 
 embryonic stem cells ( hesc ) or human induced pluripotent stem cells , offer unique possibilities of studying the very 
 early steps of human development that lead to the formation of germ 
 layers and primordial tissues .
this opportunity was seized by the 
 european commission - funded research consortium for the use of embryonic 
 stem cell - based novel alternative tests ( esnats ) for the prediction 
 of toxicity of drug candidates ( www.esnats.eu ) . within 
 this context
, several tests have been established that recapitulate 
 critical processes of early neuronal development in vitro.(4,1623 ) recently , two positive control compounds of developmental neurotoxicity , 
 valproic acid ( vpa ) and methylmercury , have been tested in four test 
 systems .
the use of two compound 
 concentrations indicated that the number of altered genes and the 
 extent of deregulation strongly depend on the concentration of the 
 test compound .
the highest concentration used in this previous study 
 was a benchmark concentration ( bmc ) that reduced overall cell viability 
 by not more than 10% .
this procedure was chosen to avoid testing of 
 cytotoxic concentrations that might generate unspecific gene expression 
 patterns due to cell death .
however , our previous study also showed 
 that exposure to a concentration that reduces viability by exactly 
 10% is difficult to perform .
however , it has 
 not yet been tested if cytotoxic concentrations really cause a death 
 associated expression signature and how the transcriptome changes 
 are affected qualitatively by concentrations deviating from the bmc .
comparisons between test systems and test compounds in this earlier 
 study suggested that different toxicants 
 may regulate the same master transcription factors .
as such factors 
 define cell identity and they show coordinated 
 regulation in different pathologies , it may be 
 assumed that they also play a role in developmental toxicity , and 
 the threshold of their activation may determine the threshold of teratogenicity .
to address the above issues , we used the ukn1 test system for more 
 detailed experiments . in this assay ,
hesc differentiate 
 to neuroepithelial precursor cells ( nep ) within 6 days , and they were 
 exposed to vpa during the entire period of differentiation .
we used 
 in this study eight concentrations of vpa covering a range from completely 
 nontoxic to strongly cytotoxic and analyzed genome wide expression 
 patterns .
vpa was chosen as a well - studied positive control compound 
 of developmental neurotoxicity .
the data were analyzed 
 by classical biostatistical approaches , but we also developed quantitative 
 measures of developmental disturbance , based on systems biology considerations . 

for instance , we quantified not only alterations of individual gene 
 ontologies ( go ) , as pioneered by the piersma lab , but we also addressed superordinate cell biological processes to 
 gain overall insight .
we report that concentrations at the uppermost 
 limit of the noncytotoxic range may be a reasonable compromise for 
 gene array classification studies because with such concentrations , 
 ( i ) the genes that indicate a teratogenic mode of action are sufficiently 
 deregulated , and ( ii ) dilution by unspecific cytotoxicity associated 
 genes only start to emerge .
gelatin , putrescine , selenium , progesterone , 
 apotransferin , glucose , insulin , and valproic acid were obtained from 
 sigma ( steinheim , germany ) .
fgf-2 ( basic fibroblast growth factor ) and noggin were obtained from 
 r&d systems ( minneapolis , mn , usa ) .
y-27632 , sb-43154 , and dorsomorphine 
 dihydrochloride were from tocris bioscience ( bristol , uk ) .
all cell culture reagents 
 were from gibco / invitrogen ( darmstadt , germany ) unless otherwise specified .
human embryonic stem 
 cells ( hesc ) ( h9 from wicells , madison , wi , usa ) were differentiated 
 according to the protocol published by chambers and colleagues with the following modifications . instead of 
 using 500 m noggin
, we used the combination of 35 m 
 noggin and 600 nm dorsomorphine together with 10 m sb-431642 
 for dual smad inhibition as described earlier .
this was used to prevent bmp and tgf signaling , and thus to achieve 
 a highly selective neuroectodermal lineage commitment . for handling 
 details ,
treatment with valproic acid ( vpa ) 
 was done with the indicated concentrations 
 from 25 m to 1000 m vpa dissolved in pbs .
dod0 medium 
 was prepared as indicated in ref ( 28 ) and supplemented with the indicated concentrations 
 of vpa .
medium supplemented with vpa was changed to dod1 , 2 , and 
 4 . in order to determine cytotoxicity ,
after the resazurin assay , 
 the medium was removed , and the cells were lysed in rna protect solution 
 ( quiagen ) .
affymetrix chip - based dna microarray analysis ( human genome 
 u133 plus 2.0 arrays ) was performed exactly as described earlier .
( 24 ) the following analyses 
 were performed 
 using the statistical programming language r - version 3.0.1 . for the 
 normalization of the entire set of 30 affymetrix gene expression arrays , 

the extrapolation strategy ( rma+ ) algorithm was used that applies background correction , log2 transformation , 
 quantile normalization , and a linear model fit to the normalized data 
 in order to obtain a value for each probe set ( ps ) on each array . 
 as reference ,
after normalization , at each concentration 
 the difference between gene expression and corresponding controls 
 was calculated ( paired design ) .
differential 
 expression was calculated using the r package limma . here , the combined information of the complete set of genes 
 is used by an empirical bayes adjustment of the variance estimates 
 of single genes .
this form of a moderated t test 
 is abbreviated here as the limma t test .
the resulting p - values were multiplicity - adjusted to control the false 
 discovery rate ( fdr ) by the benjamini  yekutieli procedure . as a result , for each concentration , a gene 
 list was obtained , with corresponding estimates for log fold change 
 and p - values of the limma t test 
 ( unadjusted and fdr - adjusted ) .
color encodes the magnitude of 
 the values , ranging from blue ( low ) to yellow ( high ) .
principal component 
 analysis ( pca ) plots were used to visualize expression data in two 
 dimensions , representing the first two principal components , that 
 is , the two orthogonal directions of the data with highest variance . 

the software r , version 3.0.1 , was used for all calculations and display 
 of pca and heat maps ( r_development_core_team 2012 ) .
the calculation and display of toxicity curves was done 
 using graphpad prism 5.0 ( graphpad software , la jolla , usa ) .
the venn 
 diagrams for the comparison of gene expression , gene ontology ( go ) 
 terms , and transcription factor binding sites ( tfbs ) between test 
 systems were constructed according to chow and rodgers .
the size of circles and areas was chosen proportional 
 to the number of elements included if possible , i.e. , if there is 
 no empty section ( zero in venn diagram ) . before clustering
first , we filtered out all of the ps with 
 too small variation . for this , we used a standard deviation cutoff 
 of 0.25 .
second , we standardized the expression measurements for each 
 ps , by subtracting its mean and dividing with its standard deviation . 

the clustering itself was performed using the k - means 
 function from r. for the index 
 of overrepresented gos , the significantly regulated genes were identified 
 for each drug concentration .
then , the sum of ogo among up - regulated 
 ps ( at a given drug concentration ) , of ogo among down - regulated ps , 
 and among all regulated ps was calculated .
this sum was divided by 
 2 ( to account for redundancies in this sum ) to arrive at the used 
 index number . for the gene regulation index ,
then , 
 for each condition , the top 100 genes were selected ( highest fold 
 changes , irrespective of the direction of regulation ) .
for each of 
 these genes , the negative log of the fdr corrected p - value was determined , and then , the sum of these 100 values was 
 formed .
the teratogenicity index is expressed by two coordinates in 
 the plane formed by the index of ogos and the gene regulation index .
the gene ontology group 
 enrichment was performed using r - version 3.0.1 with the topgo package , and only results from the biological process 
 ontology were kept .
the category was considered only if it was annotated 
 with more than 100 genes , and the minimal enrichment p - value across all concentrations was below 0.001 . a new biostatistical 
 method for measuring activation profiles of go
the method uses a segmentation 
 test to identify significant enrichment 
 with up- or down - regulated genes ( restricted by p - value < 0.05 ) , comparing each concentration with the controls . 

the large number of enrichment tests is fdr - adjusted according to 
 the two - step procedure by benjamini et al . to limit the fdr to 1% .
the profile 0++++++ means that the corresponding gene set 
 is significantly enriched with up - regulated genes at the second lowest 
 concentration and for all higher concentrations .
analogously , a gene 
 set with profile 00000 - - is enriched with down - regulated genes only 
 at the two highest concentrations .
transcription factor binding 
 site enrichment ( tfbs ) was performed 
 using the prima algorithm ( http://acgt.cs.tau.ac.il/prima/ ) provided in the expander software 
 suite ( version 6.04;http://acgt.cs.tau.ac.il / expander/ ) as described in ref ( 24 ) .
a transcription factor was included in heat maps or venn diagrams 
 only if the minimal enrichment p - value across all 
 concentrations was below 0.001 .
in order to investigate concentration dependent 
 gene expression changes induced by vpa , a recently established hesc - based in vitro system , ukn1 , was applied .
this test system recapitulates 
 the stage of human neuroectodermal induction leading to the formation 
 of neural ectodermal progenitor cells ( nep ) .
these developmental steps take place during a six days process , 
 which is covered by exposure to the test compound ( figure 1a ) . in order to define a concentration range for 
 gene array analysis ,
viability was tested in a concentration range 
 from 25 to 1000 m vpa . to this end , differentiating cells were 
 treated with the indicated concentrations of vpa , and a resazurin 
 viability assay was performed .
concentrations up to 550 m vpa 
 appeared to be noncytotoxic , whereas viability was decreased by 2030% 
 with 800 and 1000 m ( figure 1b ) .
higher 
 concentrations could not be used , as otherwise mrna extractions were 
 not feasible anymore .
on the basis of the cytotoxicity tests , these 
 eight concentrations ( 25 , 150 , 350 , 450 , 550 , 650 , 800 , and 1000 m ) 
 of vpa and untreated controls were analyzed by dna microarray ( dma ) 
 analyses ( affymetrix ) .
exposure conditions exactly the same as those 
 for the cytotoxicity experiments were used ( figure 1a ) .
embryonic stem cells ( hesc ) were converted 
 in six days of differentiation ( dod6 ) to a pure population of neuroepithelial 
 precursors ( nep ) . adapted from ref ( 24 ) .
the cells were differentiated and treated as indicated in a. on dod6 , 
 resazurin reduction was measured , and viability is given as a percentage 
 of untreated controls .
( c ) cells as treated in b were used for whole transcriptome analysis . 

the result is displayed in the form of a two - dimensional principal 
 component ( pc ) analysis diagram .
each point represents one experiment 
 (= data from one microarray ) , and the color coding indicates the concentration 
 of vpa ( in m ) used in the experiment .
all data sets ( columns ) 
 were sorted by similarity clustering . the absolute gene expression 
 levels ( log2-scaling ) of the 1000 transcripts with highest variance 
 were color - coded for display . principal component analysis ( pca ) illustrates a convincing 
 concentration 
 progression model ( figure 1c ) .
the data of 
 three independent experiments per concentration clustered closely 
 together and moved into the direction of the first principal component 
 with increasing concentrations .
an exception was the concentration 
 of 650 m vpa , which showed a high degree of variability ( figure 1c ) .
accordingly , we observed high variability in the resazurin 
 assay between the different experiments ( figure 1b ) .
therefore , this concentration was excluded from the further analysis 
 shown in the main figures , but evaluations with 650 m are available 
 in the supporting information ( table s1 ) . 

further , we generated heat maps based on z - scores and corresponding 
 dendrograms ( figure s1 , supporting information ) calculated by hierarchical clustering ( figure 1d ) .
this analysis confirmed the results from the pca ( figure 1c ) that the different experimental replicates cluster 
 closely together , apart from 650 m that was already identified 
 as an outlier . for the following analysis of differential gene 
 expression , the 
 benjamini - yekutieli false discovery rate adjustment ( by - fdr ) of the p - values
this method has the advantage that the 
 results do not vary if one concentration is added to or left out from 
 the analysis .
however , this approach can introduce a bias since by - fdr 
 is more liberal in the case of large numbers of differentially expressed 
 genes . to get an overview over the extent of the bias , the number 
 of differentially up - regulated and down - regulated genes was calculated 
 for three types of adjustment ( no adjustment , per - concentration fdr 
 adjustment , and overall fdr adjustment ) ( figure s2 , supporting information ) .
as the differences between per - concentration 
 and overall adjustment were not substantial , the analysis throughout 
 the manuscript is based on per concentration adjustment .
one reason 
 for the small differences is that genes were selected not only based 
 on p - value cut - offs but also required a > 2-fold 
 change .
after fdr correction , no significantly altered ps were 
 found for 
 the two lowest tested vpa concentrations of 25 and 150 m ( figure 2a ) . between 350 and 1000 m vpa , the numbers 
 of up- or down - regulated probe sets ( ps ) gradually increased ( figure 2a and table s1 , supporting information ) . from these data , it is not clear , whether genes regulated , e.g. , 
 at 450 m , were not at all regulated at lower concentrations 
 or whether they were regulated but did not pass our stringent selection 
 threshold ( p < 0.05 and fold change > 2 ) .
therefore , 
 we selected the ps that were regulated by 450 m vpa and determined 
 their average fold change ( fc ) . therefore , we determined the fc of 
 the same set of ps at lower and higher concentrations ( figure 2b ) .
as a complementary approach , we selected all 
 ps regulated significantly ( p < 0.05 ) at 450 m , 
 with or without a 2-fold cutoff of the fc .
then , we examined the distribution 
 of the p - values of all these ps at lower drug concentrations 
 ( figure s3 , supporting information ) . from 
 these additional statistical analyses ,
we conclude that a subgroup 
 of ps is only regulated at higher drug concentrations and not at all 
 at low exposure ; another subgroup is regulated concentration dependently , 
 i.e. , to a lower extent at low drug concentrations and to a higher 
 extent at high concentrations of vpa .
concentration - dependence of expression 
 changes and biological grouping 
 of probe sets after exposure to vpa .
( a ) number of differentially 
 expressed probe sets ( ps ) compared to the untreated control ( p < 0.05 , by - corrected ; fold change ( fc ) > 2 ) .
left 
 panel 
 ( red ) represents the amount of up - regulated ps and the right panel 
 ( green ) the amount of down - regulated ps .
( b ) the ps identified in 
 a to be regulated by 450 m vpa ( up , n = 554 ; 
 down , n = 110 ) were examined for their regulation 
 at all drug concentrations .
the average fc of this set of up - regulated 
 ps is shown in the left panel for different vpa concentrations ; the 
 data for the down - regulated ps are in the right panel .
( c ) significantly 
 up- and down - regulated ps and their potential upstream regulators 
 were identified by bioinformatic approaches .
the numbers of ogos ( left 
 panel ) and tfbs ( right panel ) are displayed for 350 , 550 , and 1000 
 m vpa .
red represents the go ( left panel ) and tfbs ( right panel ) 
 overrepresented among up - regulated ps and green from down - regulated 
 ps .
( d ) the ogos determined in c were classified in superordinate 
 cell biological processes by expert knowledge .
the identified processes 
 are given in the middle ; the left panel refers to up- and the right 
 panel down - regulated ps .
next , we analyzed the significantly ( p - value 
 < 0.05 
 after fdr correction ) regulated ps for overrepresented gene ontologies 
 terms ( ogo ) . for this , we used the differentially regulated ps of 
 350 ( low ) , 550 ( medium ) , and 1000 ( high ) m vpa .
this go enrichment 
 analysis and comparison across drug concentrations revealed that also 
 on the level of ogos the number increased with increasing vpa concentrations 
 ( figure 2c and table s2 , supporting information ) .
next , the identified ogos were further 
 assigned to superordinate cell biological processes as detailed in supporting information ( table s2 ) .
for the up - regulated 
 ps , it became clear that all identified processes were already changed 
 with low vpa concentration .
in addition , only a few more ogos were 
 identified additionally with higher concentrations ( figure 2c ) , and these additional ogos all fell into the 
 same superordinate cell biological processes ( figure 2d ) .
surprisingly , this finding applies also to stress pathways 
 and the process apoptosis / death .
for example , cell death processes 
 were overrepresented already at low concentrations , and there was 
 no new ogo emerging in this category at very high concentrations . 
 for the down - regulated ps
, few superordinate cell biological processes 
 were identified , and these processes only start to be down - regulated 
 with medium vpa concentration .
especially , neuronal 
 pathways were down - regulated from 550 m on . in summary
, these 
 data suggest that key biological processes start to be altered already 
 at low and medium drug concentrations , and hardly any further processes 
 are changed at high ( cytotoxic ) concentrations .
to investigate 
 this finding by a different approach , we generated 
 venn diagrams from the data on concentration progression .
this data 
 display illustrated that most genes altered at lower concentrations 
 are also altered at higher concentrations ( figure 3a and table s3 , supporting information ) .
more details on the concentration progression are shown in dendrograms 
 of up- and down - regulated genes ( figure s1 , supporting 
 information ) : at higher concentrations , additional genes were 
 deregulated , and many of the genes deregulated at low concentrations 
 to a small extent show lower p - values ( stronger deregulation ) 
 at higher concentrations . over 50% of the differentially expressed 
 ps were induced with high drug concentrations only ( figure 3b and table s4 , supporting information )
this situation was different for overrepresented transcription 
 factor binding sites ( tfbs ) ( figure 3b and 
 table s4 , supporting information ) .
a relatively 
 large fraction of all tfbs found in this study at any concentration 
 was already overrepresented at low drug levels of 350 m .
this 
 was particularly obvious for the overrepresented tfbs of down - regulated 
 genes , illustrated in figure 3b .
in summary , 
 the regulation of ps and tfbs showed a characteristically different 
 behavior : while many ps were regulated at high vpa concentrations 
 only , the corresponding tfbs were already changed at low concentrations .
( a ) the number of significantly up ( upper panel)- 
 and down - regulated ( lower panel ) ps ( cut off : p - value 
 < 0.05 , by corrected and fc > 2 ) for each concentration were 
 compared 
 to the ps of the next higher concentration .
( b ) comparison of up- 
 and down - regulated ps ( and of tfbs overrepresented among these ps ) 
 between the indicated concentrations .
overlaps are displayed as venn 
 diagrams with absolute numbers of altered elements in the circles , 
 and the number of nonaltered elements in the lower right corner .
the 
 percentage number refers to elements affected by 1000 m only , 
 relative to all affected elements .
also , ogo between different vpa concentrations were compared 
 in 
 a venn diagram to visualize the degree of overlap .
this analysis showed 
 that ogo groups behave more similarly to tfbs than to ps ( figure 4 and table s5 , supporting information ) : for up- and down - regulated genes , most ogo groups started to be 
 overrepresented at 350 m , and relatively few ogo groups are 
 added at higher concentrations .
these findings are illustrated at 
 higher detail in dendrograms ( figure s1 , supporting 
 information ) .
a relatively large fraction of ogo groups overrepresented 
 at 350 m overlaps with ogo groups overrepresented at 550 and 
 1000 m ( n = 55 ) . for down - regulated genes , 
 ogo groups began to emerge at slightly higher concentrations , with 
 only 6 ogo groups being overrepresented at 350 m .
the overall 
 features remained the same : 22 of the 25 go overrepresented at 1000 
 m were already overrepresented at 550 m ( figure 4b ) .
again , we also classified the ogos in superordinate 
 cell biological processes . also , on these higher levels , the results 
 were confirmed .
taken together , these analyses illustrate different 
 progression models for genes on the one side and ogo groups and tfbs 
 on the other . increasing vpa concentrations
however , most ogo groups 
 and tfbs were already overrepresented at 350 or 450 m with 
 relatively little added at higher concentrations .
( a ) up - regulated and ( b ) down - regulated ps were analyzed for ogo , 
 and then the identified ogo were compared between drug concentrations . 

the overlaps of ogo at vpa concentrations of 350 m ( low ) , 550 
 m ( medium ) , and 1000 m are displayed as venn diagrams . 

right : the ogos from specified areas of the diagrams ( see legend ; 
 550  1000 means overlap area of 550 m and 1000 m 
 circles ) were grouped according to superordinate cell biological processes , 
 and the absolute numbers within these groups are displayed .
the previous analysis focused on the differentially 
 expressed ps for each concentration compared to the control .
, we grouped the responses into clusters 
 of 5 different up - regulation responses and 5 types of down - regulations 
 ( figure 5a and table s6 , supporting information ) .
the number of ps falling into each 
 cluster and the number of ogos among these ps are given in figure 5b .
response courses showed 
 a monotonic pattern , but clusters 1 , 5 , and 6 contained ps with nonmonotonic 
 regulation patterns .
together , the latter clusters comprised 6.8% 
 of all ps ( n = 9624 ) considered in the cluster analysis .
( b ) the number of ps in each 
 cluster ( solid bars ) and ogos within each cluster ( hatched bars ) are 
 displayed .
red , up - regulated ; green , down - regulated ; gray , nonmonotonic . 
 ( c , d ) all ogos containing less than 1000 genes are displayed for 
 clusters 5 and 10 .
clusters 5 and 10 were 
 chosen for further investigation , as they 
 contained the genes that were regulated at high drug concentrations 
 ( 800 and 1000 m ) only .
nevertheless , a relatively high number 
 of go terms was overrepresented among these genes ( figure 5b and table s7 , supporting information ) .
the ogo identified were relatively heterogeneous , and no dominant 
 superordinate cell biological process was identified . in particular , 
 there was no overrepresentation of genes relating to cell death , cell 
 stress , or degeneration ( figure 5c , d ) .
this 
 result was well in line with our earlier data . to further follow up 
 on this , and to investigate which types of responses were triggered 
 at which drug concentration , we chose in the following analyses a 
 more targeted approach to specifically investigate the concentration 
 dependency of the regulation of biological processes .
all analyses in 
 the previous paragraphs ( figure 2 to figure 5 ) are driven by single gene results .
first , genes 
 were combined to clusters with respect to similar expression trajectories 
 across the whole concentration range
. then functional gene groups 
 ( e.g. , go groups ) were analyzed for overrepresentation with genes 
 from these identified clusters .
a meaningful alternative is to use 
 a go group ( or any other functional gene group representing , e.g. , 
 a biological process ) as the starting point for the analysis .
therefore , 
 we applied a new biostatistical method that analyzes the activity 
 ( up- or down - regulation ) of genes annotated to a go group across all 
 concentrations .
the method uses a segmentation test to identify significant enrichment with up- or down - regulated 
 genes ( restricted by p - value < 0.05 ) , comparing 
 each concentration with the controls .
for example , the profile 
 0++++++ means that the corresponding gene set is significantly enriched 
 with up - regulated genes at the second lowest concentration and for 
 all higher concentrations but not at the lowest concentration .
analogously , 
 a gene set with profile 00000 - - is enriched with down - regulated genes 
 only at the two highest concentrations but not at the five lower concentrations . 

the result of all go activation profiles is given in supporting information ( table s8 ) . for
visualization of the 
 go activation profiles , we calculated a quantitative go activation 
 score and plotted the concentration  response profiles for some 
 exemplary gos ( figure 6a ) .
next , we classified 
 the ogos identified with this method in superordinate biological processes 
 as described before .
then , transcriptome analysis was 
 performed , and differentially expressed genes were identified as in 
 figure 1 .
the data sets were examined for overrepresentation 
 of any go , and go - based concentration 
response profiles were 
 generated as detailed in table s8 ( supporting 
 information ) .
( a ) quantitative go activation score 
 was calculated by multiplying the percentage of genes within the go 
 that was found to be significantly regulated with the average fold 
 change of these regulations .
positive values reflect ogo among up - regulated 
 genes , and negative values indicate ogo among down - regulated genes . 

note the 
 different concentrations at which a go is first turned on and the 
 different types of concentration  response behavior at high 
 concentrations .
all go profiles among up- and down - regulated 
 genes were sorted , and the following three groups were selected for 
 further analysis : low ( go found to be up - regulated at all vpa concentrations 
 starting from 150 m , clusters 0++++++ and 00+++++ ) ; intermediate 
 ( go found to be up - regulated at all vpa concentrations higher than 
 350 m , clusters 000++++ and 0000+++ ) ; and high ( go found to 
 be up - regulated only at 800 and 1000 m , clusters 00000++ and 
 000000 + ) .
then , the go were assigned to superordinate cell biological 
 processes ( e.g. , migration and adhesion ) , and we calculated how many 
 of all ogo in the low  medium 
a color code was applied 
 with no color for values < 10% , yellow for 1020% , orange 
 for 2030% , and red for > 30% .
the concentration dependence of the regulation of superordinate 
 biological processes indicated three major groups : those mainly affected 
 at high concentrations ( 00000++ and 000000 + ) , those activated already 
 at medium concentrations and then being maintained ( 000++++ and 0000+++ ) , 
 and finally those represented already among low concentrations and 
 remaining up - regulated at higher concentrations ( 0++++++ and 00+++++ ) 
 ( figure 6b ) .
these findings agreed well with 
 those obtained by other analytical approaches , such as unbiased k - means 
 clustering ( figure 5 ) or comparison of the 
 concentration dependence of ogo and ps ( figures 2 and 4 ) .
application of the go profiler analysis 
 to the concentration data illustrates that only two key biological 
 processes ( catabolism and cell division ) were overrepresented exclusively 
 at cytotoxic concentrations ( 800 and 1000 m ) without being 
 enriched in the noncytotoxic range ( figure 6b and table s8 , supporting information ) .
go groups involved in catabolism of dna and proteins were overrepresented 
 for the up- and down - regulated genes at cytotoxic concentrations . 

moreover , cell division associated genes were overrepresented for 
 the down - regulated genes only at 800 and 1000 m vpa .
the superordinate 
 biological processes representing best the known mechanism of action 
 of vpa , protein acetylation and epigenetic processes , were found ( as 
 expected ) to start at noncytotoxic concentrations ( figure 6b and table s8 , supporting information ) .
another process clearly activated at noncytotoxic concentrations 
 is migration and adhesion ( figure 6b and table 
 s8 , supporting information ) .
other biological 
 processes that are of interest for the mechanism of action of vpa 
 are development and differentiation , metabolism , signal transduction , 
 and response to external compound .
go groups belonging to these processes 
 were overrepresented in the noncytotoxic but also the cytotoxic concentration 
 range .
response 
 behavior of the regulation of genes vs go and superordinate cell biological 
 processes suggested that these features may be used to develop a quantitative 
 measure of teratogenic activity .
for this purpose , the effect of each 
 drug concentration was defined by two measures , the gene regulation 
 index and the go overrepresentation index . when these values were 
 plotted onto a coordinate system formed by the two indices
, it became 
 evident that the distance from the origin was a useful measure of 
 teratogenicity ( figure 7a ) .
this way of data 
 presentation also clearly revealed that vpa concentrations between 
 125 and 450 m resulted in a progressive deviation of normal 
 differentiation ( increasing numbers of ogo among deregulated ps ) . 
 at higher concentrations , the extent of gene deregulation increased , 
 but ogos remained fairly constant .
the concentrations marking the 
 distinct increase of teratogenicity according to this presentation 
 correlate well with the vpa plasma concentrations associated with 
 human or animal birth defects . the curve formed by the teratogenicity 
 index data ( figure 7a ) differed fundamentally 
 from the cytotoxicity curve ( figure 1b ) .
the 
 latter showed strong and progressive changes at high drug concentrations , 
 while the former did not change at high vpa levels but rather at medium 
 levels .
thus , it seems feasible to define a specific teratogenicity 
 measure that yields clearly different information from plain cytotoxicity 
 in the same assay .
( a ) transcriptome data were obtained for multiple 
 drug ( vpa ) concentrations as in figure 1 . to 
 obtain a measure of the deviation of the transcriptome from normal 
 ( teratogenicity measure ) ,
two indices were calculated and used to 
 define the two dimensions of a developmental toxicity plane .
the index is the 
 sum of the negative logarithms of the p - values of 
 the 100 most regulated genes .
the second dimension reflects the extent 
 of coordinated changes in biological processes reflected by gos .
the 
 index is proportional to the number of overrepresented gos in the 
 gene sets .
( b ) summary of overall findings on concentration - dependent 
 transcriptome deviations : drug concentrations were chosen in a way 
 to allow either normal neuroectodermal differentiation of human embryonic 
 stem cells ( hesc ) or disturbed differentiation ( teratogenic concentration 
 range ) or cytotoxicity . over this large concentration range ( 251000 
 m ) , the number of deregulated genes increased continuously , 
 once a certain threshold concentration was reached ( 125 m ) . 

in contrast to this , superordinate biological regulations , as indicated 
 by enriched gos ( gene ontologies ) or tfbs ( transcription factor binding 
 sites ) , increased steeply in the teratogenic range and then more or 
 less reached a plateau . at cytotoxic concentrations ,
thus , cytotoxicity , overall gene 
 regulation , and coordinate regulation of biological processes showed 
 largely different concentration dependencies .
there were clear lower 
 thresholds , and the extent of transcriptome regulation as such , without 
 further bioinformatic analysis , did not appear to be a good measure 
 for teratogenicity .
in 
 this study , we used a stem cell based test system that recapitulates 
 neural differentiation of hesc . in order to elucidate the optimal 
 concentration for toxicogenomics
, we challenged this test system with 
 increasing concentration of the well - known neurodevelopmental toxicants 
 vpa .
we found that the number of differentially regulated ps continuously 
 increased but that up - targets like tfbs and functional gene ontologies 
 reach saturation already at the lower border of intermediate concentrations . 

in addition , we show that with clearly cytotoxic concentrations only 
 few more superordinate biological processes occur .
therefore , we conclude 
 that the most sensitive concentration for dnt toxicogenomic testing 
 would be the highest nontoxic concentration of a compound
. however , 
 small deviations from this concentration will not disrupt the results 
 on the basis of ogo and tfbs .
we present here the development of the 
 transcriptomics - based teratogenicity index that exactly reflects this 
 behavior for the comparison of deregulated ps and ogos ( figure 7a ) .
although transcriptomics is frequently applied 
 in toxicology , it has remained unclear how it is best used for toxicological 
 statements and predictions .
basic issues , such as the question of 
 appropriate concentrations , timing of exposure , and timing of sampling , 
 have not been addressed in a systematic and quantitative manner .
only 
 very few studies , mainly in the field of carcinogenesis , have attempted 
 to correlate the extent or type of overall transcriptome changes with 
 toxicological predictions in a quantitative manner .
for instance , 
 thomas and colleagues have correlated benchmark doses 
 based on go terms and classical pathology with cancer development . 

we are not aware of any such study in developmental toxicology , and 
 our study is a first attempt in this direction , even though only a 
 single model compound , vpa , has been studied . given the astonishing dearth of data in the field , 
 we feel that our study is an important beginning but should not be 
 overinterpreted .
further data are required to see whether the response 
 features observed here are representative of a larger group of developmental 
 toxicants .
not every cell system is equally suited for studies 
 on design principles 
 of gene regulation .
our test system provided 
 the required homogeneity and synchronicity to perform such analysis . 

previously others and we have shown that the differentiation protocol 
 delivers a > 90% cell population of the same cell type . moreover , genome wide mrna analysis demonstrated a monotonic progression 
 model of deregulated 
 genes .
a very high fraction of genes deregulated at a lower concentration 
 were also deregulated at the next higher concentration .
moreover , 
 the next higher concentration usually yielded additional genes that 
 have not yet been deregulated at the lower concentration .
these are 
 good conditions to study the control mechanisms and principles of 
 the cellular response to increasing concentrations .
unsupervised 
 clustering indicated that 
 there is a smaller subfraction of ps that behaved differently .
to 
 follow up on this , we selected the ps that were regulated by an intermediate 
 drug concentration ( 450 m ) and determined their average fold 
 change ( fc ) . then , we determined the fc of the same set of ps at lower 
 and higher concentrations . at the lowest concentration ( 25 m )
the linear middle part supports the suggestion that the selected 
 ps were regulated increasingly more with increasing concentration , 
 and this may have affected the number that was detected after filtering 
 for significance and fc .
the flattening of the curve on the lower 
 and upper ends suggests that most likely additional mechanisms are 
 involved , i.e. , regulation of different genes at different concentrations , 
 as observed here by other methods . as a complementary approach , we 
 selected the ps that were significantly regulated at 450 m 
 ( p < 0.05 ) and calculated the p - values for the same ps at low concentrations .
most of the ps lost 
 significance at 150 m and nearly all at 25 m .
the situation 
 was somewhat different , when the initial selection of ps included 
 a fc cutoff of 2 .
under this condition , all ps selected at 450 m 
 still reached a p < 0.05 significance level at 
 350 m , and most of them still did so at 150 m . at 25 
 m , still about 10% of the ps were significant .
this suggests , 
 that some ps that we considered here as regulated only at high drug 
 concentrations but not at low drug exposure were indeed already affected 
 at low concentrations but that the change was too small to be picked 
 up as significant .
this feature may have to be considered for future 
 systems biology based quantitative models of gene regulation .
comprehensive biostatistical analyses including the establishment 
 of a novel gene ontology activation profiler ( goap ) suggest the following 
 concept of three concentration ranges .
a range of tolerance is observed up to 25 m vpa . at 25 m ,
not a single 
 up- or down - regulated gene with p < 0.05 and a 
 fold change higher than two was obtained .
this observation is in agreement 
 with the general concept that toxicity based on , e.g. , enzyme inhibition 
 or induction , or receptor activation , acts by threshold mechanisms . in the cytotoxic concentration range ( 800 m 
 and 1000 m )
, we expected that an apoptosis signature may occur , 
 but we observed the induction of apoptosis , death , and stress associated 
 ogo groups already at the intermediate concentration range . also , 
 go analysis of the genes that are exclusively up- or down - regulated 
 with the cytotoxic concentrations revealed no additional death related 
 ogos .
in addition , goap analysis revealed that go that are related 
 to catabolism and metabolism are induced with cytotoxic concentrations . 
 among the down - regulated ps , cell division
analysis of gos overrepresented among ps in this range gave evidence 
 for developmental disturbances , cell migration , and the down - regulation 
 of neuronal pathways .
therefore , we defined 
 the range from 150 m to 550 m vpa as the deregulated / teratogenic concentration range ( figure 7b ) .
these findings 
 are in good agreement with human data that indicate that vpa plasma 
 concentrations above 500 m are associated with strongly increasing 
 risk for malformations .
the newly established 
 goap method revealed additional and more sensitive answers for this 
 range of teratogenic concentrations . between 150 and 550 m , 
 a number of go groups became overrepresented that reflect the specific 
 mode of action of vpa , namely , protein acetylation and epigenetic 
 processes .
down - regulated ps comprises genes involved in all kind 
 of epigenetic processes like several 
 histone acetyl transferases ( e.g. , kat2b , kat6a , kat5 ) , polycomb proteins 
 ( e.g. , suz12 ) , lysine methyl transferases ( e.g. , ehmt2 ) , and hdacs 
 themselves ( e.g. , hdac 2 , 5 , and 9 ) .
regarding go regulation 
 patterns over the concentration range , 
 we also identified go groups that showed , e.g. , a 00+++00 nonmonotonic 
 pattern .
. we could verify the main 
 monotonic regulation pattern by another analysis targeting upstream 
 regulators of transcriptional changes .
we found that most tfbs are already induced 
 with teratogenic concentrations , and only few additional tfbs were 
 induced with cytotoxic concentrations ( figure 7b ) .
these comprised ap-2 and pax5 for up - regulated and ipf1 and pbx 
 for down - regulated ps .
on the basis of differentially expressed 
 ps and altered functional 
 processes , we developed a teratogenicity index that could distinguish 
 between dnt - relevant ( the teratogenic range ) and cytotoxic concentrations . 

for this purpose , the extent of transcriptome deregulation was captured 
 by measures for deregulation of individual genes and a measure for 
 biological processes ( ogo ) .
although cytotoxicity only starts at 800 m 
 vpa , the functional processes were significantly changed already with 
 350 m and the teratogenicity index only slightly changed with 
 800 m and 1000 m . therefore , we present here an index 
 that distinguishes between noncytotoxic , teratogenic , and cytotoxic 
 concentrations and could help in the future to improve the experimental 
 design of toxicogenomic studies . in particular , our data support the 
 use of the highest noncytotoxic concentration of toxicogenomics studies , 
 if no other information on the test compounds is available . at this 
 concentration ,
a maximum power is reached ( number of deregulated ps ) 
 without compromising the biological information obtained ( e.g. , altered 
 ogo or tfbs ) .
accidental use of slightly too high concentrations would 
 be expected to have only minor effects on the overall outcome . in conclusion , the data from this concentration progression study 
 suggest a steady increase of deregulated genes .
however , they mostly 
 belong to the same biological processes and seem to be controlled 
 by the same tf that are already deregulated at the intermediate concentrations 
 of 350 and 450 m vpa ( figure 7b ) .
superenhancers are unusual enhancers , which are found at genes responsible 
 for cell type identity and occupied by key transcription factors for 
 this cell type . in esc ,
e.g. , these dna sequences are mainly occupied 
 by the key stemness factors oct4 , nanog , and sox2 .
it was also shown that these superenhancers play crucial 
 roles in the development of many diseases .
therefore , it is very likely that they are also disturbed or differentially 
 occupied by a toxic compound and play a crucial role in the interpretation 
 of our upstream target analysis . in the future , it will be interesting 
 to investigate whether deviations from the normal development on the 
 level of tfbs and superordinate cell biological processes will be 
 shared by diverse toxicants and whether the time of toxicant exposure 
 has an effect on such findings .
 information on design principles 
 governing transcriptome changes 
 upon transition from safe to hazardous drug concentrations or from 
 tolerated to cytotoxic drug levels are important for the application 
 of toxicogenomics data in developmental toxicology . here 
 , we tested 
 the effect of eight concentrations of valproic acid ( vpa ; 251000 
 m ) in an assay that recapitulates the development of human 
 embryonic stem cells to neuroectoderm . 
 cells were exposed to the drug 
 during the entire differentiation process , and the number of differentially 
 regulated genes increased continuously over the concentration range 
 from zero to about 3000 . 
 we identified overrepresented transcription 
 factor binding sites ( tfbs ) as well as superordinate cell biological 
 processes , and we developed a gene ontology ( go ) activation profiler , 
 as well as a two - dimensional teratogenicity index . 
 analysis of the 
 transcriptome data set by the above biostatistical and systems biology 
 approaches yielded the following insights : ( i ) tolerated ( 25 
 m ) , deregulated / teratogenic ( 150550 m ) , and 
 cytotoxic ( 800 m ) concentrations could be differentiated . 
 
 ( ii ) biological signatures related to the mode of action of vpa , such 
 as protein acetylation , developmental changes , and cell migration , 
 emerged from the teratogenic concentrations range . ( 
 iii ) cytotoxicity 
 was not accompanied by signatures of newly emerging canonical cell 
 death / stress indicators , but by catabolism and decreased expression 
 of cell cycle associated genes . 
 ( iv ) most , but not all of the go groups 
 and tfbs seen at the highest concentrations were already overrepresented 
 at 350450 m . 
 ( v ) the teratogenicity index reflected 
 this behavior , and thus differed strongly from cytotoxicity . 
 our findings 
 suggest the use of the highest noncytotoxic drug concentration for 
 gene array toxicogenomics studies , as higher concentrations possibly 
 yield wrong information on the mode of action , and lower drug levels 
 result in decreased gene expression changes and thus a reduced power 
 of the study . 
Introduction
Materials and Methods
Results
Discussion