Regional variation in neurovascular coupling and why we still lack a Rosetta Stone

Functional magnetic resonance imaging (fMRI) is the dominant tool in cognitive neuroscience although its relation to underlying neural activity, particularly in the human brain, remains largely unknown. A major research goal, therefore, has been to uncover a ‘Rosetta Stone’ providing direct translation between the blood oxygen level-dependent (BOLD) signal, the local field potential and single-neuron activity. Here, I evaluate the proposal that BOLD signal changes equate to changes in gamma-band activity, which in turn may partially relate to the spiking activity of neurons. While there is some support for this idea in sensory cortices, findings in deeper brain structures like the hippocampus instead suggest both regional and frequency-wise differences. Relatedly, I consider four important factors in linking fMRI to neural activity: interpretation of correlations between these signals, regional variability in local vasculature, distributed neural coding schemes and varying fMRI signal quality. Novel analytic fMRI techniques, such as multivariate pattern analysis (MVPA), employ the distributed patterns of voxels across a brain region to make inferences about information content rather than whether a small number of voxels go up or down relative to baseline in response to a stimulus. Although unlikely to provide a Rosetta Stone, MVPA, therefore, may represent one possible means forward for better linking BOLD signal changes to the information coded by underlying neural activity.This article is part of the theme issue ‘Key relationships between non-invasive functional neuroimaging and the underlying neuronal activity’.     

Negative BOLD fMRI response in the visual cortex carries precise stimulus-specific information

Sustained positive BOLD (blood oxygen level-dependent) activity is employed extensively in functional magnetic resonance imaging (fMRI) studies as evidence for task or stimulus-specific neural responses. However, the presence of sustained negative BOLD activity (i.e., sustained responses that are lower than the fixation baseline) has remained more difficult to interpret. Some studies suggest that it results from local "blood stealing" wherein blood is diverted to neurally active regions without a concomitant change of neural activity in the negative BOLD regions. However, other evidence suggests that negative BOLD is a result of local neural suppression. In both cases, regions of negative BOLD response are usually interpreted as carrying relatively little, if any, stimulus-specific information (hence the predominant reliance on positive BOLD activity in fMRI). Here we show that the negative BOLD response resulting from visual stimulation can carry high information content that is stimulus-specific. Using a general linear model (GLM), we contrasted standard flickering stimuli to a fixation baseline and found regions of the visual cortex that displayed a sustained negative BOLD response, consistent with several previous studies. Within these negative BOLD regions, we compared patterns of fMRI activity generated by flickering Gabors that were systematically shifted in position. As the Gabors were shifted further from each other, the correlation in the spatial pattern of activity across a population of voxels (such as the population of V1 voxels that displayed a negative BOLD response) decreased significantly. Despite the fact that the BOLD signal was significantly negative (lower than fixation baseline), these regions were able to discriminate objects separated by less than 0.5 deg (at approximately 10 deg eccentricity). The results suggest that meaningful, stimulus-specific processing occurs even in regions that display a strong negative BOLD response.     

Confidence Intervals for fMRI Activation Maps

Neuroimaging activation maps typically color voxels to indicate whether the blood oxygen level-dependent (BOLD) signals measured among two or more experimental conditions differ significantly at that location. This data presentation, however, omits information critical for interpretation of experimental results. First, no information is represented about trends at voxels that do not pass the statistical test. Second, no information is given about the range of probable effect sizes at voxels that do pass the statistical test. This leads to a fundamental error in interpreting activation maps by naïve viewers, where it is assumed that colored, “active” voxels are reliably different from uncolored “inactive” voxels. In other domains, confidence intervals have been added to data graphics to reduce such errors. Here, we first document the prevalence of the fundamental error of interpretation, and then present a method for solving it by depicting confidence intervals in fMRI activation maps. Presenting images where the bounds of confidence intervals at each voxel are coded as color allows readers to visually test for differences between “active” and “inactive” voxels, and permits for more proper interpretation of neuroimaging data. Our specific graphical methods are intended as initial proposals to spur broader discussion of how to present confidence intervals for fMRI data.     

Measuring Granger Causality between cortical regions from voxelwise fMRI BOLD signals with LASSO

Functional brain network studies using the Blood Oxygen-Level Dependent (BOLD) signal from functional Magnetic Resonance Imaging (fMRI) are becoming increasingly prevalent in research on the neural basis of human cognition. An important problem in functional brain network analysis is to understand directed functional interactions between brain regions during cognitive performance. This problem has important implications for understanding top-down influences from frontal and parietal control regions to visual occipital cortex in visuospatial attention, the goal motivating the present study. A common approach to measuring directed functional interactions between two brain regions is to first create nodal signals by averaging the BOLD signals of all the voxels in each region, and to then measure directed functional interactions between the nodal signals. Another approach, that avoids averaging, is to measure directed functional interactions between all pairwise combinations of voxels in the two regions. Here we employ an alternative approach that avoids the drawbacks of both averaging and pairwise voxel measures. In this approach, we first use the Least Absolute Shrinkage Selection Operator (LASSO) to pre-select voxels for analysis, then compute a Multivariate Vector AutoRegressive (MVAR) model from the time series of the selected voxels, and finally compute summary Granger Causality (GC) statistics from the model to represent directed interregional interactions. We demonstrate the effectiveness of this approach on both simulated and empirical fMRI data. We also show that averaging regional BOLD activity to create a nodal signal may lead to biased GC estimation of directed interregional interactions. The approach presented here makes it feasible to compute GC between brain regions without the need for averaging. Our results suggest that in the analysis of functional brain networks, careful consideration must be given to the way that network nodes and edges are defined because those definitions may have important implications for the validity of the analysis.     

Validity and Calibration of the Youth Activity Profile

Purpose

The purpose of this study was to calibrate and cross-validate the Youth Activity Profile (YAP), a self-report tool designed to capture physical activity (PA) and sedentary behaviors (SB) in youth.

Methods

Eight schools in the Midwest part of the U.S. were involved and a total of 291 participants from grades 4–12 agreed to wear an accelerometer (SWA Armband) and complete the YAP in two separate weeks (5–7 days apart). Individual YAP items capture PA behavior during specific segments of the week and these items were combined with temporally matched estimates of moderate-to-vigorous PA (MVPA) and sedentary time from the SWA to enable calibration. Quantile regression procedures yielded YAP prediction algorithms that estimated MVPA at School, MVPA at Out-of-School, MVPA on Weekend, as well as time spent in SB. The YAP estimates of time spent in MVPA and SB were cross-validated using Pearson product correlations and limits of agreement, as indicative of individual error and, equivalence testing techniques as indicative of group-level error.

Result

Following calibration, the correlations between YAP and SWA estimates of MVPA were low to moderate (r_range = .19 to .58) and individual-level YAP estimates of MVPA ranged from -134.9% to +110.0% of SWA MVPA values. Differences between aggregated YAP and SWA MVPA ranged from -3.4 to 21.7 minutes of MVPA at the group-level and predicted YAP MVPA estimates were within 15%, 20%, and 30%, of values from the SWA for the School, Out-of-School, and Weekend time periods, respectively. Estimates of time spent in SB were highly correlated with each other (r = .75). The individual estimates of SB ranged from -54.0% to +44.0% of SWA sedentary time, and the aggregated group-level estimates differed by 49.7 minutes (within 10% of the SWA aggregated estimates).

Conclusions

This study provides preliminary evidence that the calibration procedures enabled the YAP to provide estimates of MVPA and SB that approximated values from an objective monitor. The YAP provides a simple, low-cost and educationally sound method to accurately estimate children’s MVPA and SB at the group level.     

Decoding Multiple Sound Categories in the Human Temporal Cortex Using High Resolution fMRI

Perception of sound categories is an important aspect of auditory perception. The extent to which the brain’s representation of sound categories is encoded in specialized subregions or distributed across the auditory cortex remains unclear. Recent studies using multivariate pattern analysis (MVPA) of brain activations have provided important insights into how the brain decodes perceptual information. In the large existing literature on brain decoding using MVPA methods, relatively few studies have been conducted on multi-class categorization in the auditory domain. Here, we investigated the representation and processing of auditory categories within the human temporal cortex using high resolution fMRI and MVPA methods. More importantly, we considered decoding multiple sound categories simultaneously through multi-class support vector machine-recursive feature elimination (MSVM-RFE) as our MVPA tool. Results show that for all classifications the model MSVM-RFE was able to learn the functional relation between the multiple sound categories and the corresponding evoked spatial patterns and classify the unlabeled sound-evoked patterns significantly above chance. This indicates the feasibility of decoding multiple sound categories not only within but across subjects. However, the across-subject variation affects classification performance more than the within-subject variation, as the across-subject analysis has significantly lower classification accuracies. Sound category-selective brain maps were identified based on multi-class classification and revealed distributed patterns of brain activity in the superior temporal gyrus and the middle temporal gyrus. This is in accordance with previous studies, indicating that information in the spatially distributed patterns may reflect a more abstract perceptual level of representation of sound categories. Further, we show that the across-subject classification performance can be significantly improved by averaging the fMRI images over items, because the irrelevant variations between different items of the same sound category are reduced and in turn the proportion of signals relevant to sound categorization increases.     

Large upward bias in estimation of locus-specific effects from genomewide scans

The primary goal of a genomewide scan is to estimate the genomic locations of genes influencing a trait of interest. It is sometimes said that a secondary goal is to estimate the phenotypic effects of each identified locus. Here, it is shown that these two objectives cannot be met reliably by use of a single data set of a currently realistic size. Simulation and analytical results, based on variance-components linkage analysis as an example, demonstrate that estimates of locus-specific effect size at genomewide LOD score peaks tend to be grossly inflated and can even be virtually independent of the true effect size, even for studies on large samples when the true effect size is small. However, the bias diminishes asymptotically. The explanation for the bias is that the LOD score is a function of the locus-specific effect-size estimate, such that there is a high correlation between the observed statistical significance and the effect-size estimate. When the LOD score is maximized over the many pointwise tests being conducted throughout the genome, the locus-specific effect-size estimate is therefore effectively maximized as well. We argue that attempts at bias correction give unsatisfactory results, and that pointwise estimation in an independent data set may be the only way of obtaining reliable estimates of locus-specific effect-and then only if one does not condition on statistical significance being obtained. We further show that the same factors causing this bias are responsible for frequent failures to replicate initial claims of linkage or association for complex traits, even when the initial localization is, in fact, correct. The findings of this study have wide-ranging implications, as they apply to all statistical methods of gene localization. It is hoped that, by keeping this bias in mind, we will more realistically interpret and extrapolate from the results of genomewide scans.     

From Blood Oxygenation Level Dependent (BOLD) Signals to Brain Temperature Maps

A theoretical framework is presented for converting Blood Oxygenation Level Dependent (BOLD) images to brain temperature maps, based on the idea that disproportional local changes in cerebral blood flow (CBF) as compared with cerebral metabolic rate of oxygen consumption (CMRO ₂) during functional brain activity, lead to both brain temperature changes and the BOLD effect. Using an oxygen limitation model and a BOLD signal model, we obtain a transcendental equation relating CBF and CMRO ₂ changes with the corresponding BOLD signal, which is solved in terms of the Lambert W function. Inserting this result in the dynamic bioheat equation describing the rate of temperature changes in the brain, we obtain a nonautonomous ordinary differential equation that depends on the BOLD response, which is solved numerically for each brain voxel. Temperature maps obtained from a real BOLD dataset registered in an attention to visual motion experiment were calculated, obtaining temperature variations in the range: (−0.15, 0.1) which is consistent with experimental results. The statistical analysis revealed that significant temperature activations have a similar distribution pattern than BOLD activations. An interesting difference was the activation of the precuneus in temperature maps, a region involved in visuospatial processing, an effect that was not observed on BOLD maps. Furthermore, temperature maps were more localized to gray matter regions than the original BOLD maps, showing less activated voxels in white matter and cerebrospinal fluid.     

Functional Connectivity in MRI Is Driven by Spontaneous BOLD Events

Functional brain signals are frequently decomposed into a relatively small set of large scale, distributed cortical networks that are associated with different cognitive functions. It is generally assumed that the connectivity of these networks is static in time and constant over the whole network, although there is increasing evidence that this view is too simplistic. This work proposes novel techniques to investigate the contribution of spontaneous BOLD events to the temporal dynamics of functional connectivity as assessed by ultra-high field functional magnetic resonance imaging (fMRI). The results show that: 1) spontaneous events in recognised brain networks contribute significantly to network connectivity estimates; 2) these spontaneous events do not necessarily involve whole networks or nodes, but clusters of voxels which act in concert, forming transiently synchronising sub-networks and 3) a task can significantly alter the number of localised spontaneous events that are detected within a single network. These findings support the notion that spontaneous events are the main driver of the large scale networks that are commonly detected by seed-based correlation and ICA. Furthermore, we found that large scale networks are manifestations of smaller, transiently synchronising sub-networks acting dynamically in concert, corresponding to spontaneous events, and which do not necessarily involve all voxels within the network nodes oscillating in unison.     

Dynamic Adjustment of Stimuli in Real Time Functional Magnetic Resonance Imaging

The conventional fMRI image analysis approach to associating stimuli to brain activation is performed by carrying out a massive number of parallel univariate regression analyses. fMRI blood-oxygen-level dependent (BOLD) signal, the basis of these analyses, is known for its low signal-noise-ratio and high spatial and temporal signal correlation. In order to ensure accurate localization of brain activity, stimulus administration in an fMRI session is often lengthy and repetitive. Real-time fMRI BOLD signal analysis is carried out as the signal is observed. This method allows for dynamic, real-time adjustment of stimuli through sequential experimental designs. We have developed a voxel-wise sequential probability ratio test (SPRT) approach for dynamically determining localization, as well as decision rules for stopping stimulus administration. SPRT methods and general linear model (GLM) approaches are combined to identify brain regions that are activated by specific elements of stimuli. Stimulus administration is dynamically stopped when sufficient statistical evidence is collected to determine activation status across regions of interest, following predetermined statistical error thresholds. Simulation experiments and an example based on real fMRI data show that scan volumes can be substantially reduced when compared with pre-determined, fixed designs while achieving similar or better accuracy in detecting activated voxels. Moreover, the proposed approach is also able to accurately detect differentially activated areas, and other comparisons between task-related GLM parameters that can be formulated in a hypothesis-testing framework. Finally, we give a demonstration of SPRT being employed in conjunction with a halving algorithm to dynamically adjust stimuli.     

Pairwise statistical significance of local sequence alignment using sequence-specific and position-specific substitution matrices

Pairwise sequence alignment is a central problem in bioinformatics, which forms the basis of various other applications. Two related sequences are expected to have a high alignment score, but relatedness is usually judged by statistical significance rather than by alignment score. Recently, it was shown that pairwise statistical significance gives promising results as an alternative to database statistical significance for getting individual significance estimates of pairwise alignment scores. The improvement was mainly attributed to making the statistical significance estimation process more sequence-specific and database-independent. In this paper, we use sequence-specific and position-specific substitution matrices to derive the estimates of pairwise statistical significance, which is expected to use more sequence-specific information in estimating pairwise statistical significance. Experiments on a benchmark database with sequence-specific substitution matrices at different levels of sequence-specific contribution were conducted, and results confirm that using sequence-specific substitution matrices for estimating pairwise statistical significance is significantly better than using a standard matrix like BLOSUM62, and than database statistical significance estimates reported by popular database search programs like BLAST, PSI-BLAST (without pretrained PSSMs), and SSEARCH on a benchmark database, but with pretrained PSSMs, PSI-BLAST results are significantly better. Further, using position-specific substitution matrices for estimating pairwise statistical significance gives significantly better results even than PSI-BLAST using pretrained PSSMs.     

A New Test for Selection Applied to Codon Usage in Drosophila simulans and D. mauritiana

In many organisms, synonymous codon usage is biased by a history of natural selection. However, codon bias, itself, does not indicate that selection is ongoing; it may be a vestige of past selection. Simple statistical tests have been devised to infer ongoing selection on codon usage by comparing the derived state frequency spectra at polymorphic sites segregating either derived preferred codons or derived unpreferred codons; if selection is effective, the frequency of derived states should be higher in the former. We propose a new test that uses the inferred degree of preference, essentially calculating the correlation of derived state frequency and the difference in preference between the derived and the ancestral states; the correlation should be positive if selection is effective. When implementing the test, derived and ancestral states can be assigned by parsimony or on the basis of relative probability. In either case, statistical significance is estimated by a simple permutation test. We explored the statistical power of the test by sampling polymorphism data from 14 loci in 16 strains of D. simulans, finding that the test retains 80% power even when quite a few of the data are discarded. The power of the test likely reflects better use of multiple features of the data, combining population frequencies of polymorphic variants and quantitative estimates of codon preferences. We also applied this novel test to 14 newly sequenced loci in five strains of D. mauritiana, showing for the first time ongoing selection on codon usage in this species.     

Reproducibility of swallow-induced cortical BOLD positive and negative fMRI activity

Functional MRI (fMRI) studies have demonstrated that a number of brain regions (cingulate, insula, prefrontal, and sensory/motor cortices) display blood oxygen level-dependent (BOLD) positive activity during swallow. Negative BOLD activations and reproducibility of these activations have not been systematically studied. The aim of our study was to investigate the reproducibility of swallow-related cortical positive and negative BOLD activity across different fMRI sessions. We studied 16 healthy volunteers utilizing an fMRI event-related analysis. Individual analysis using a general linear model was used to remove undesirable signal changes correlated with motion, white matter, and cerebrospinal fluid. The group analysis used a mixed-effects multilevel model to identify active cortical regions. The volume and magnitude of a BOLD signal within each cluster was compared between the two study sessions. All subjects showed significant clustered BOLD activity within the known areas of cortical swallowing network across both sessions. The cross-correlation coefficient of percent fMRI signal change and the number of activated voxels across both positive and negative BOLD networks were similar between the two studies (r ≥ 0.87, P < 0.0001). Swallow-associated negative BOLD activity was comparable to the well-defined "default-mode" network, and positive BOLD activity had noticeable overlap with the previously described "task-positive" network. Swallow activates two parallel cortical networks. These include a positive and a negative BOLD network, respectively, correlated and anticorrelated with swallow stimulus. Group cortical activity maps, as well as extent and amplitude of activity induced by volitional swallowing in the cortical swallowing network, are reproducible between study sessions.     

Taxon sampling, correlated evolution, and independent contrasts

Independent contrasts are widely used to incorporate phylogenetic information into studies of continuous traits, particularly analyses of evolutionary trait correlations, but the effects of taxon sampling on these analyses have received little attention. In this paper, simulations were used to investigate the effects of taxon sampling patterns and alternative branch length assignments on the statistical performance of correlation coefficients and sign tests; "full-tree" analyses based on contrasts at all nodes and "paired-comparisons" based only on contrasts of terminal taxon pairs were also compared. The simulations showed that random samples, with respect to the traits under consideration, provide statistically robust estimates of trait correlations. However, exact significance tests are highly dependent on appropriate branch length information; equal branch lengths maintain lower Type I error than alternative topological approaches, and adjusted critical values of the independent contrast correlation coefficient are provided for use with equal branch lengths. Nonrandom samples, with respect to univariate or bivariate trait distributions, introduce discrepancies between interspecific and phylogenetically structured analyses and bias estimates of underlying evolutionary correlations. Examples of nonrandom sampling processes may include community assembly processes, convergent evolution under local adaptive pressures, selection of a nonrandom sample of species from a habitat or life-history group, or investigator bias. Correlation analyses based on species pairs comparisons, while ignoring deeper relationships, entail significant loss of statistical power and as a result provide a conservative test of trait associations. Paired comparisons in which species differ by a large amount in one trait, a method introduced in comparative plant ecology, have appropriate Type I error rates and high statistical power, but do not correctly estimate the magnitude of trait correlations. Sign tests, based on full-tree or paired-comparison approaches, are highly reliable across a wide range of sampling scenarios, in terms of Type I error rates, but have very low power. These results provide guidance for selecting species and applying comparative methods to optimize the performance of statistical tests of trait associations.     

Echo-time and field strength dependence of BOLD reactivity in veins and parenchyma using flow-normalized hypercapnic manipulation

While the BOLD (Blood Oxygenation Level Dependent) contrast mechanism has demonstrated excellent sensitivity to neuronal activation, its specificity with regards to differentiating vascular and parenchymal responses has been an area of ongoing concern. By inducing a global increase in Cerebral Blood Flow (CBF), we examined the effect of magnetic field strength and echo-time (TE) on the gradient-echo BOLD response in areas of cortical gray matter and in resolvable veins. In order to define a quantitative index of BOLD reactivity, we measured the percent BOLD response per unit fractional change in global gray matter CBF induced by inhaling carbon dioxide (CO(2)). By normalizing the BOLD response to the underlying CBF change and determining the BOLD response as a function of TE, we calculated the change in R(2)(*) (ΔR(2)(*)) per unit fractional flow change; the Flow Relaxation Coefficient, (FRC) for 3T and 1.5T in parenchymal and large vein compartments. The FRC in parenchymal voxels was 1.76±0.54 fold higher at 3T than at 1.5T and was 2.96±0.66 and 3.12±0.76 fold higher for veins than parenchyma at 1.5T and 3T respectively, showing a quantitative measure of the increase in specificity to parenchymal sources at 3T compared to 1.5T. Additionally, the results allow optimization of the TE to prioritize either maximum parenchymal BOLD response or maximum parenchymal specificity. Parenchymal signals peaked at TE values of 62.0±11.5 ms and 41.5±7.5 ms for 1.5T and 3T, respectively, while the response in the major veins peaked at shorter TE values; 41.0±6.9 ms and 21.5±1.0 ms for 1.5T and 3T. These experiments showed that at 3T, the BOLD CNR in parenchymal voxels exceeded that of 1.5T by a factor of 1.9±0.4 at the optimal TE for each field.     

BOLD Granger Causality Reflects Vascular Anatomy

A number of studies have tried to exploit subtle phase differences in BOLD time series to resolve the order of sequential activation of brain regions, or more generally the ability of signal in one region to predict subsequent signal in another region. More recently, such lag-based measures have been applied to investigate directed functional connectivity, although this application has been controversial. We attempted to use large publicly available datasets (FCON 1000, ADHD 200, Human Connectome Project) to determine whether consistent spatial patterns of Granger Causality are observed in typical fMRI data. For BOLD datasets from 1,240 typically developing subjects ages 7–40, we measured Granger causality between time series for every pair of 7,266 spherical ROIs covering the gray matter and 264 seed ROIs at hubs of the brain’s functional network architecture. Granger causality estimates were strongly reproducible for connections in a test and replication sample (n=620 subjects for each group), as well as in data from a single subject scanned repeatedly, both during resting and passive video viewing. The same effect was even stronger in high temporal resolution fMRI data from the Human Connectome Project, and was observed independently in data collected during performance of 7 task paradigms. The spatial distribution of Granger causality reflected vascular anatomy with a progression from Granger causality sources, in Circle of Willis arterial inflow distributions, to sinks, near large venous vascular structures such as dural venous sinuses and at the periphery of the brain. Attempts to resolve BOLD phase differences with Granger causality should consider the possibility of reproducible vascular confounds, a problem that is independent of the known regional variability of the hemodynamic response.     

Inferring species richness using multispecies occupancy modeling: Estimation performance and interpretation

Multispecies occupancy models can estimate species richness from spatially replicated multispecies detection/non‐detection survey data, while accounting for imperfect detection. A model extension using data augmentation allows inferring the total number of species in the community, including those completely missed by sampling (i.e., not detected in any survey, at any site). Here we investigate the robustness of these estimates. We review key model assumptions and test performance via simulations, under a range of scenarios of species characteristics and sampling regimes, exploring sensitivity to the Bayesian priors used for model fitting. We run tests when assumptions are perfectly met and when violated. We apply the model to a real dataset and contrast estimates obtained with and without predictors, and for different subsets of data. We find that, even with model assumptions perfectly met, estimation of the total number of species can be poor in scenarios where many species are missed (>15%–20%) and that commonly used priors can accentuate overestimation. Our tests show that estimation can often be robust to violations of assumptions about the statistical distributions describing variation of occupancy and detectability among species, but lower‐tail deviations can result in large biases. We obtain substantially different estimates from alternative analyses of our real dataset, with results suggesting that missing relevant predictors in the model can result in richness underestimation. In summary, estimates of total richness are sensitive to model structure and often uncertain. Appropriate selection of priors, testing of assumptions, and model refinement are all important to enhance estimator performance. Yet, these do not guarantee accurate estimation, particularly when many species remain undetected. While statistical models can provide useful insights, expectations about accuracy in this challenging prediction task should be realistic. Where knowledge about species numbers is considered truly critical for management or policy, survey effort should ideally be such that the chances of missing species altogether are low.     

Selection of the experimental unit in teratology studies.

In teratology experiments the litter (pregnant female) rather than the fetus is advocated as being the proper experimental unit upon which to base the statistical analysis. It is pointed out that per litter tests, by being based on the average fetal response within a litter, do take the individual fetus into account. Actual experimental data are used to show that when litter effects are present a per fetus analysis is invalid and may seriously exaggerate the significance level. It is also shown that there appears to be little loss in sensitivity in performing per litter tests even in the unlikely event that there are no litter effects. Thus it certainly seems prudent to analyze teratology data with test procedures that treat the litter as the experimental unit.