### Biological experimental design

Primary human epidermal keratinocytes (HEK, Cascade Biologics) were cultured with serum-free defined media at 5% CO_{2} and 37°C. The HEK were exposed to the following low-micron to nano-scale materials: carbonyl iron (FC), carbon black (CB) silica (SiO_{2}) and single-walled carbon nanotubes (SWNT). FC was supplied by ISP Technologies with a mean particle size of 5.8 μm, CB was supplied by Degussa with a mean particle size of 17 nm, and SiO_{2} was supplied by US Silica with a mean particle size of 1.6 μm. SWNT was prepared by SouthWest Nanotechnologies with a mean diameter of 0.8 nm. The purity of SWNT was checked, and it was found that all heavy metal contamination was very low, less than 1%. The nanotubes were acid-purified from the raw manufactured intermediate. The average length of the SWNT nanotubes was 960 nm; this is the only nanomaterial among the ones considered here that is not spherical. From previous viability assays [13, 15], two different doses were selected for the exposure of each nanomaterial – a noncytotoxic dose and a cytotoxic dose. The cytotoxic dose was extrapolated from the viability curve of each substance where 50% of the cells were still viable. The noncytotoxic and cytotoxic doses for each nanomaterial were found to be, respectively, 0.03 and 1 (FC), 0.01 and 0.5 (CB), 0.1 and 1 (SiO_{2}), and 0.001 and 1 (SWNT), all concentrations being expressed in mg/ml. The cells were exposed to FC, CB, SiO_{2}, and SWNT and harvested at 2, 4, 6, 8, 12, 18, and 24 hours after exposure. Cells cultured under the same conditions and not exposed to any of the nanomaterials were harvested at the same time points for a time-matched control baseline.

The cells were trysinized and cell counts taken. The cells were then collected by centrifugation, snap frozen in liquid nitrogen, and stored at -80°C. The frozen pellets were shipped to GenUs Biosystems for isolation of the total RNA and processing of the gene expression microarrays. Total RNA was isolated from the cell pellets, reverse-transcribed to biotinylated cRNAs and hybridized onto 10 K human gene expression microarrays (GE Healthcare). The corresponding cRNA for each biological sample was hybridized to triplicate microarrays. The arrays were rinsed, dried, scanned and image analyzed and the flat files returned to the Houston Advanced Research Center.

### Data analysis and preprocessing

The microarray flat files contained the quantitative expression values for all probes (positive controls, negative controls, fiducial, and discovery) and all discovery probe values were assessed with a quality flag. Only those probes that have "DISCOVERY" and "G" (Good-quality) tags across all time points and dose levels of all treatments were considered for further analysis to guarantee maximal reproducibility of the results. By using this criterion, the number of genes was reduced from 10,458 to 2,370. Then the average of the three replicate measurements was considered as the actual expression level of each gene at that time in its corresponding treatment.

### Heatmaps

The heatmaps are visualizations of hierarchical clustering [16] using "average" linkage and Pearson Correlation as the distance metric between expression levels of each gene after taking the log ratio of the data. It is worthy to mention that Pearson correlation is in fact the normalized correlation between values of two random variables that have had their mean subtracted from them. The dendrogram added to the left side of each heat map is obtained by Pearson correlation. The heatmaps were drawn by first subtracting the mean of each row from each data cell and then normalizing each row to obtain a variance of 1. Note that normalization was performed after taking the log ratio of the expression of each gene between the given treatment and the control baseline. Heatmaps obtained in this way reveal the underlying shape of expression pattern of genes better and are less affected by some large values which may be the result of noise. For more information on applying hierarchical clustering to gene expression patterns, see [16] and [17] and the references mentioned there.

### Multi-dimensional Scaling (MDS)

Multidimensional Scaling (MDS) is a method to project high dimensional data onto a lower dimension while maintaining the approximate distances between data points [21]. The accuracy of the representation is measured by the "stress." A stress of 10% is usually considered good.

First, the data are collected by taking the log of the expression values of genes in each treatment and the baseline. Then at each time point, a matrix of distances between the gene-expression profiles corresponding to each treatment and the baseline is created, where the distance is calculated as 1 – the Pearson correlation between the profiles. Finally, MDS is performed based on this matrix, in order to obtain a low-dimensional representation of the gene-expression profiles and the distance among them.

### Self-Organizing Maps (SOM)

SOM [15] is a powerful method for clustering. After log-normalizing the data as explained above, this method was employed in two ways.

First, the underlying hidden number of clusters of the data was estimated using two different criteria the silhouette criterion [18] and the CLEST criterion [19]. In [19], it has been shown that the CLEST criterion is the stronger tool to find number of clusters. However, from the comparison of the CLEST and the silhouette criteria in [19], it is inferred that whenever the true number of clusters has been equal to 2, the silhouette criterion is the same as the CLEST or has done a better job than the CLEST criterion, while for greater number of clusters, the CLEST performs better. Therefore, the criterion that we have used to determine the number of clusters is a combination of the two criteria. First, the silhouette criterion is implemented to determine the number of clusters. If the number of clusters for any treatment is equal to 2, and if the silhouette width for this number of clusters is fairly close to 1 (e.g. greater than 0.6) then we assigned 2 as the optimal number of clusters. Otherwise, the CLEST criterion is implemented to determine the number of clusters.

In the second way of implementing SOM, larger numbers of clusters than the previous approach were tried, as a means of obtaining clusters that show tighter expression patterns than in the first approach.

### ANOVA model

A two-way ANOVA model was used to identify genes for which there are significant differences in the response to the various dose levels, and have them ranked by the *p*-values of dose main effect. The general procedure was as follows:

0 – For each gene, compile the data in a two-dimensional array consisting of the expression values from the three levels of dose (untreated, low and high) and the eight time points. Use Tukey's additivity test for two-way ANOVA (one observation per cell) to see if there is interaction between time and dose for the given gene. If there is no interaction go to step 2.

1 – If there is significant interaction between these two factors, try to remove the interaction by using the Box-Cox transformation [20]. The removal of interaction can be tested by Tukey's additivity test for transformed values corresponding to each parameter of the Box-Cox transformation. We have chosen this parameter to vary between -2 to 2 with 0.1 spacing. If the interaction is removed, use these as the new data for the given gene. If the interaction could not be removed, set the p-value of dose main effect to NA and go to step 0 to start the process for the next gene.

2 – Use two-way ANOVA to find whether there are significant differences among the levels of the dose and time factor. If there is not any significant difference on any factor, go to step 4.

3 – Use Tukey's HSD post-hoc test to determine between which levels of either factor there are significant differences and record them. Since Tukey's HSD test is conservative, in rare cases the test cannot find any significant difference between levels.

4 – Save the p-value for the difference in dose main effect for this gene and go to step 0 to start the process for the next gene.

At the end of this process, sort all genes with respect to p-values on dose main effect. The lower the p-value is for a gene, the stronger the evidence is that there are significant differences in response among the various levels of dose for that gene.