Depth relationships and measures of tissue thickness in dorsal midbrain.

Dorsal human midbrain contains two nuclei with clear laminar organization, the superior and inferior colliculi. These nuclei extend in depth between the superficial dorsal surface of midbrain and a deep midbrain nucleus, the periaqueductal gray matter (PAG). The PAG, in turn, surrounds the cerebral aqueduct (CA). This study examined the use of two depth metrics to characterize depth and thickness relationships within dorsal midbrain using the superficial surface of midbrain and CA as references. The first utilized nearest-neighbor Euclidean distance from one reference surface, while the second used a level-set approach that combines signed distances from both reference surfaces. Both depth methods provided similar functional depth profiles generated by saccadic eye movements in a functional MRI task, confirming their efficacy for delineating depth for superficial functional activity. Next, the boundaries of the PAG were estimated using Euclidean distance together with elliptical fitting, indicating that the PAG can be readily characterized by a smooth surface surrounding PAG. Finally, we used the level-set approach to measure tissue depth between the superficial surface and the PAG, thus characterizing the variable thickness of the colliculi. Overall, this study demonstrates depth-mapping schemes for human midbrain that enables accurate segmentation of the PAG and consistent depth and thickness estimates of the superior and inferior colliculi.

INTRODUCTION

Human midbrain mediates a panoply of critical brain functions that range from homeostasis (Edlow, McNab, Witzel, & Kinney, 2016; Hommel et al., 2006; Tsujino & Sakurai, 2009) to perception (Behbehani, 1995; Sprague, 1991) to cognition (Morgane, Galler, & Mokler, 2005; Roeper, 2013). However, studying the human dorsal midbrain, including the superior colliculus (SC) and inferior colliculus (IC), is challenging due to the small size and deep location of its component nuclei, as well as limited methods for understanding their internal structure. As a result, most of our existing knowledge about the function of the midbrain arises from research conducted in animal models. As human ultra‐high field MRI adoption grows, the constraints posed by collicular size and location are reduced due to increased MR signal at 7 T and above. However, the challenge of understanding the internal structure of the midbrain remains, requiring a depth‐mapping scheme for probing the anatomy and function of the human midbrain.

Both SC and IC have an intricate laminar cytoarchitecture and exist in a complex 3D topology that is folded into four small hillocks on the dorsal surface of the midbrain. The inner boundaries of the colliculi abut another nucleus, the periaqueductal gray (PAG), which surrounds the ventricular cerebral aqueduct (CA).

The laminar functional organization of mammalian colliculi has been clearly established based on several electrophysiology and lesion studies in animal models. Experiments in nonhuman primates and cats have shown that the layers of SC can be divided into three groups with independent functional purposes. Superficial layers have been shown to receive direct visual input and have retinotopically organized receptive fields (Cynader & Berman, 1972; Feldon & Kruger, 1970). Intermediate layers mediate oculomotor control (Robinson, 1972), while deep layers are associated with multimodal inputs, primarily multisensory, and visuomotor neurons (Meredith & Stein, 1986; Sprague & Meikle, 1965). Likewise, electrophysiological studies in IC of rats and primates have shown that the central nucleus organizes auditory inputs of varying frequency in a laminar fashion (Baumann et al., 2010; Schreiner & Langner, 1997). Invasive studies have exhibited this tonotopic organization along the dorsal‐to‐ventral direction, which roughly corresponds to laminar depth (Baumann et al., 2011; Cheung et al., 2012; Loftus, Malmierca, Bishop, & Oliver, 2008; Malmierca et al., 2008; Schreiner & Langner, 1997).

While functional MRI (fMRI) has been used extensively in human cerebral cortex, the human brainstem has been relatively neglected despite its critical role in brain function. Functional contrast‐to‐noise ratio (CNR) is often low due to its deep location (Singh, Pfeuffer, Zhao, & Ress, 2017) and functional data is susceptible to partial volume effects due to the small size of important brainstem structures compared with conventional fMRI voxel sizes of 3–4 mm. However, several recent high‐resolution fMRI studies using millimeter‐scale voxels have demonstrated detailed analysis of the functional laminar topography of the human colliculi with reasonable CNR. Indeed, these studies have confirmed the functional properties of the three major layer groups of SC (Katyal & Ress, 2014; Linzenbold & Himmelbach, 2012; Loureiro et al., 2017; Savjani, Katyal, Halfen, Kim, & Ress, 2018; Schneider & Kastner, 2010), as well as the auditory organization of IC (De Martino et al., 2013; Moerel, De Martino, Uğurbil, Yacoub, & Formisano, 2015; Ress & Chandrasekaran, 2013), encouraging interest in resolving these laminar variations of human colliculi in more detailed experiments. A recent fMRI study showed depth‐dependent variations of the BOLD response with 1‐mm sampling in SC (Loureiro et al., 2017). This work utilized binary mask erosion techniques to generate depth‐dependent regions‐of‐interest (ROI). The approach resolves depth at the 1‐mm scale but does not permit associations between laminae at different depths.

In previous work, we utilized a nearest‐neighbor Euclidean point‐to‐surface definition of depth to define layers in SC by creating depth kernels that associate a particular locus of superficial tissue with deeper tissue. These kernels were small computational cylinders extending from the surface of the brainstem inwards, normal to the surface of SC. With this method, we successfully demonstrated depth‐dependent BOLD responses from visual stimulation and attention (Katyal & Ress, 2014; Katyal, Zughni, Greene, & Ress, 2010) as well as saccadic activity in SC (Savjani et al., 2018) and frequency‐dependent auditory stimulation in IC (Ress & Chandrasekaran, 2013).

However, defining a laminar neighborhood of tissue using nearest‐neighbor associations is a drastic oversimplification. In the small convoluted colliculi, the cylinders are susceptible to oversampling or undersampling deep tissue, depending on the gradient of the curvature. In order to capture the physical structure of collicular layers as observed in histology and MR microscopy, a more topologically consistent method of defining depth within the colliculi is desirable. Several methods have been developed to compute depth in cortex, including solutions of Laplace's equation (Jones, Buchbinder, & Aharon, 2000) and equi‐volume methods (Waehnert et al., 2014). However, they are computationally expensive to solve in a stable and accurate fashion, and have never been applied to midbrain.

Here, we compare two methods for depth‐analysis in human midbrain: our previous Euclidean approach, and an algebraic level‐set algorithm that utilizes two surfaces. We adapted the level‐set method originally implemented in cortex (Kim, Taylor, & Ress, 2017) to human midbrain using the superficial brainstem tissue‐CSF boundary as the outer surface and the CA as the inner surface to create a depth coordinate normalized to the thickness of the tissue under investigation. First, we compared the ability of both methods to delineate visual‐ and saccade‐evoked function in superficial SC, concluding that both methods perform similarly. We also delineated the inner boundaries of the colliculi in order to set a precedent for future functional studies in the deepest layers of SC. This was accomplished by localizing the outer boundary of the PAG using a combination of Euclidean and level‐set methods applied to gray‐matter tissue‐probability maps obtained at 9.4T. Lastly, we evaluated anatomical thickness of the colliculi, using both methods to measure the distance from the collicular surface to the outer boundary of the PAG. In this case, the level‐set approach provided more precise measures of tissue depth in the superior and inferior colliculus. Altogether, we show that our level‐set method for depth analysis meets or exceeds the performance of more typical Euclidean approaches and suggests new applications of our method for mapping the human brainstem.

METHODS

Subjects

We applied our level‐set scheme to 20 subjects. Eight of these participated in subcortical fMRI experiments at 3T (Savjani et al., 2018), giving informed consent under procedures reviewed and authorized by the Baylor College of Medicine Institutional Review Board. The other 12 subjects participated in quantitative MRI experiments at 9.4T (Hagberg et al., 2017) and underwent a physical and psychological check‐up by a local physician and provided written informed consent following local research ethics policies and procedures. These investigations were conducted in agreement with the World Medical Association Declaration of Helsinki in its most recent version (2013).

Acquisition and segmentation of anatomical images

For each of the eight subjects scanned at 3T, we acquired high‐resolution (0.7‐mm voxels) T1‐weighted anatomical volumes using an MP‐RAGE sequence (TI = 900 ms, TR = 2,600 ms, 9° flip angle, TA = 22 min) on a 3T Siemens (Siemens Medical Solutions, Erlangen, Germany) Trio scanner with a product 32‐channel head coil. The remaining 12 subjects were scanned at 9.4T using a Siemens Magnetom scanner equipped with a 16‐channel, dual‐row transmit array operating in the circularly polarized mode, and a 31‐channel receive array (Shajan et al., 2014). We acquired B1 field maps by the actual flip angle imaging method (Yarnykh, 2010; Yarnykh & Yuan, 2004) with nominal flip angle FA = 60°; repetition time TR1/TR2 = 20/100 ms; echo time TE = 7 ms, voxel size = 3 × 3 × 5 mm3; and acquisition time TA = 225 s. We further acquired high‐resolution (0.8‐mm isotropic voxels) quantitative whole brain T1‐maps using an MP2‐RAGE sequence (TI1/TI2 = 900/3500 ms, TR = 6 ms, volume TR = 9 s, TE = 2.3 ms, GRAPPA = 3, partial‐Fourier factor 6/8, 256 RF pulses, and TA = 9 min 40 s). At 9.4T, this approach permitted high‐quality anatomical imaging with the much shorter acquisition time than at 3T. The MP2RAGE images were reconstructed off‐line while correcting for deviations from the nominal excitation flip angle and for T2‐dependent deviations in inversion efficiency of the adiabatic inversion pulse (Hagberg et al., 2017). From the quantitative T1 maps, two sets of synthetic, B1‐artifact‐free, T1‐weighted, MP2RAGE contrast images were generated pixel‐wise from the analytical model equations described in the Appendix of Marques et al. (2010). The first image set was used for brain tissue segmentation in volume space using SPM12. The white‐matter gray matter tissue contrast could be increased with respect to the image acquisition by setting the TI2 to 2,650 ms, and the inversion efficiency to 0.9 while the remaining parameters were the same as in the MR‐acquisition. The second image set was generated to improve the quality of the FreeSurfer segmentation (performed using 100 iterations and the following recon_all options: “–highres ‐3T –schwartzya3t‐atlas”) and was based on the following parameters: TI1/TI2 = 950/2100 ms; flip angle = 5/3°; volume TR = 6 s; and an inversion efficiency of 0.8451 (corresponding to the experimentally determined median value found in brain tissue).

For all subjects, to generate the outer segmentation, the brainstem tissue (BS) of each subject was initially identified using a probabilistic Bayesian approach in the FreeSurfer 6.0 software package (Iglesias et al., 2015). These segmentations were then edited manually to obtain precise delineation of the four colliculi and their vicinity. To create the inner segmentation, the CA was identified using a mixture of manual and automatic region‐growing tools implemented in ITK‐SNAP (Yushkevich et al., 2006).

Surface modeling

A surface model S 1 was estimated at the interface between the brainstem tissue and adjacent CSF or thalamic tissues (Figure 1a). An initial isosurface was created directly from the segmentation using MATLAB R2016a (Mathworks Corp., Natick, MA). Then, to reduce voxelation artifacts, we applied five iterations of refinement using a variational, volume‐preserving deformable surface algorithm (Bajaj, Xu, & Zhang, 2008; Khan et al., 2011; Xu, Pan, & Bajaj, 2006). The same procedure was applied to the CA segmentation to construct a second surface, S 2 (Figure 1e).

FIGURE 1

Isosurfaces of tissue proximal to the cerebral aqueduct constructed at five values of w for one subject scanned at 3T, with w = 0 representing the interface between brainstem tissue and CSF (S 1) and w = 1 representing the cerebral aqueduct (S 2)

3D level‐set depth‐mapping

We calculated a normalized signed distance function w using two separate physical distances relative to the corresponding surfaces (Figure 2a,b). First, Euclidean distances were calculated for voxels in the volume to the nearest triangle of the designated surface S and S, giving a 3D mapping of each voxel to its distance value (Eberly, 1999). Second, the sign for each voxel was determined based on its location; voxels enclosed within the surface were assigned positive distance values, while voxels outside the surface were given negative distances. Thus, distances from S 1 (d 1) became increasingly positive from the surface to BS into deeper tissue, while distances from S 2 (d 2) were only positive within CA and negative in BS and surrounding CSF. Finally, we calculated a weighted sum of the signed distance functions, with distance parameter w, using a level‐set scheme that solves an Eikonal equation (Bajaj et al., 2008; Khan et al., 2011; Kim et al., 2017):where the level‐set parameter w provides a normalized depth coordinate between the two surfaces. The depth metric w (Figure 2c) is zero on the brainstem surface () and unity at the surface of CA (); w < 0 is outside BS and w > 1 within CA. Because the normalized depth coordinate is independent of the variable tissue thickness, the above equation also enables generation of isosurfaces that evolve smoothly from S 1 to S 2 (Figure 1a–e) as well as outside of the domain bounded by the surfaces. Taking advantage of this computational feature, the changes in curvature that distinguish the colliculi became less obvious with increasing w (moving from brainstem surface to CA), and the level‐set surfaces of w ultimately morph into the cerebral aqueduct, S 2. Normalized depth maps were consistent between subjects (Figure 3), showing a smooth progression between the brainstem surface at w = 0 to the CA at w = 1. As evident in the surfaces shown in Figure 1, the characteristic curvature of the colliculi diminished with increasing w depth.

FIGURE 2

Distance maps generated for one subject overlaid on an axial T1 MRI image of superior colliculus (far left) acquired at 3T. Shown in panels to the right are Euclidean signed‐distance functions (a) d 1, (b) d 2, (c) normalized depth, w and (d) level‐set depth derived from w by ray tracing. Black dashed lines show the surface representations of the superior colliculus (S 1) in (a) and the cerebral aqueduct (S 2) in (b), at signed‐distance function values of zero

FIGURE 3

Left to right: axial, coronal, and sagittal cross sections of normalized depth maps in four subjects scanned at 3T, and location of cross sections (pink: axial, green: coronal, and blue: sagittal) on brainstem surface models (S 1) of each subject. Axial cross‐sections show both superior colliculus (Subjects 1 and 2) and inferior colliculus (Subjects 3 and 4). Sagittal cross sections ran through the cerebral aqueduct (Subject 1), through the crown of the colliculi (Subjects 2 and 3), and through tissue adjacent to the cerebral aqueduct (Subject 4)

Generation of streamlines and depth‐averaging kernels

To obtain unique correspondences between S 1 and S 2, we treated w as a pseudopotential. We calculated by convolution with five‐point stencil kernels along each dimension. Then, we used ray tracing of to generate streamlines originating at each vertex of S 1 (Figure 4a) and propagating throughout midbrain and vicinity. The streamlines were initialized along the corresponding surface normal of S 1 and propagated in small piecewise increments (0.25 voxels), following from w = 0 to w = 1.5. The same process was repeated following the negative gradient from w = 0 to w = −1. The positive and negative streamline segments were then concatenated into a continuous streamline from w = −1 to w = 1.5. To avoid numerical pathologies, we utilized several heuristics. Topological inconsistencies were avoided by terminating streamlines that encountered large changes (>80°) in the direction of in a single step. Additionally, stagnating streamlines that failed to make sufficient spatial progress after each iteration beyond a minimum threshold (typically 0.05 voxels) before reaching CA were removed. Any streamlines that remained incomplete (w < 1) were also discarded after a maximum number of iterations (typically 64 in the direction of positive w and 32 in the negative direction).

FIGURE 4

(a) Quasi‐axial view of a surface representation of superior colliculus (S 1) and cerebral aqueduct (S 2) in one subject scanned at 3T with three level‐set depth‐averaging kernels (radius = 0.7 mm, manifold distance). Their component streamlines (dark green) are regridded (light green) to the spatial resolution of the anatomical volume (0.7 mm isotropic voxels). (b) In comparison, cylindrical depth‐averaging kernels demonstrate undesirable overlap (“contention”, red arrow) and inappropriately directed sampling in deeper collicular tissue. (c) Quantification of the deviation metric, Δd, between level‐set and Euclidean sampling methods. (d) Relationship between Euclidean nearest‐neighbor depth and level‐set depth in collicular tissue. (e) Surface representations of brainstem (S 1) of one subject, with deviation Δd values at three level‐set depths mapped back to the surface vertices from which each streamline originates

To create depth‐averaging kernels that associated voxels on the surface of BS with deeper tissue, we gathered all streamlines within a chosen radius, typically 0.7 mm, of manifold distance around each vertex on S 1 (Figure 4a). The coordinates of each streamline were then rounded to the nearest voxel and associated with their corresponding vertex of origin on the brainstem tissue surface.

For comparison, we utilized our earlier approach that used only the Euclidean distance metric d 1 to define depth. In this method, depth kernels were computed by extending a 0.7‐mm manifold radius disk of BS surface voxels along their mean normal in both directions, yielding a cylinder of voxels (Figure 4b).

Depth profiles of visual stimulation‐saccade‐evoked activity in SC

To test the new level‐set streamline approach, we re‐analyzed existing data on the polar‐angle representation of saccadic eye movements in SC (Savjani et al., 2018). For both stimulation and saccades, we measured the peak BOLD response at depths between 0.5 mm outside BS and 3.5 mm inside BS at increments of 0.1 mm, averaging the values within each 1.2‐mm‐wide bin. Given that the nature of our stimulus was expected to elicit narrow bands of saccade‐evoked activity in intermediate SC, as detailed in (Savjani et al., 2018), we automatically generated elliptical ROIs at the surface of SC to capture the most robust activity. We expected the profiles for saccades to be shifted deeper into SC compared with those from visual stimulation. The reliability of our data was established using a bootstrapping method. We resampled across runs with replacement over 2000 iterations, calculating a new laminar profile for each resampled average. This data were used to obtain the mean depth profiles across subjects with their corresponding 68% confidence intervals, and depth values of peak activity with their 68% confidence intervals and p‐values (Figure 5). Depth profiles generated using streamlines and level‐set physical tissue depth were compared to profiles using our previous cylinders and Euclidean depth.

FIGURE 5

Bootstrapped (n = 2000) mean laminar profiles of normalized BOLD response in superior colliculus (SC) across 4 subjects scanned at 3T, depth‐averaged using: (a) Euclidean depth and cylinder depth kernels, (b) level‐set depth and streamline depth kernels, and (c) normalized w depth and streamline depth kernels. Shaded regions around curves and peak locations indicate 68% confidence intervals

PAG analysis in midbrain

We estimated the outer boundary of the PAG using quantitative structural data obtained at 9.4T. For this purpose, the unified segmentation approach (Ashburner & Friston, 2005) and a gray‐matter tissue‐probability atlas (Lorio et al., 2016) in SPM12 were applied to the MP2RAGE‐based image sets with optimized white‐gray matter contrast (as described above). This map was qualitatively discriminative of the PAG boundary at very high probability thresholds (Figure 6a). To quantify an optimal probability threshold (), we minimized spatial uncertainty (Ress, Harlow, Marshall, & McMahan, 2004) of the PAG boundary across subjects. Since the outer boundary of the deep PAG roughly forms an annulus surrounding CA, we chose to define it in units of our Euclidean distance metric from CA, d 2. We used our new level‐set depth kernels to quantify P as a function of depth, binning increasing values of d 2 from S 1 to S 2, at increments of 0.05 mm (Figure 6b). We initially calculated depth profiles for comparison using normalized w distance, level‐set path distance, and Euclidean distance from S 1, as well as using the cylinder depth kernels. However, they provided almost no contrast between the PAG and the surrounding tissue, and using d 2 as a reference depth metric clearly proved to be the most useful method. From the streamline depth profiles, we then computed the spatial uncertainty U:where is the SD of P across subjects at each value of d . Examining the mean values of U with respect to gray matter probability, we determined that the optimal = 0.96. We then constructed and characterized a 3D estimate of the PAG boundary (Figure 6c,d). The segmentation for CA was upsampled by a factor of two and skeletonized to create an axis. Using this axis as a normal vector, 2D slices were taken through the brainstem probability map at each point along the axis. These contours were then fit in a least‐squared error sense with ellipses centered around the cerebral aqueduct. These ellipses were characterized by three parameters: left–right axis length a, dorsal‐ventral axis length b, and aspect ratio . (Figure 6e). These 2D ellipses were then flood filled and transformed back to the locations of the slices they were calculated from, to create a 3D segmentation of the elliptical fit of the PAG. We termed the resulting surface generated from this segmentation S .

FIGURE 6

Gray‐matter probability map of one subject scanned at 9.4T with contour lines drawn at p = .96 (red). (b) Left: Individual depth profiles (gray) sampling gray‐matter probability maps using streamlines, with respect to distance from cerebral aqueduct (CA). Mean depth profiles across subjects is shown in black. Right: Individual values of spatial uncertainty U with respect to gray‐matter probability. Mean across subjects is shown in black. (c) Slices of the gray‐matter probability map of one subject taken normal to CA. Three sections are shown running through superior colliculus (SC; 1), intercollicular sulcus between SC and inferior colliculus (IC; 2), and IC (3). Contour lines drawn at p = .96 are shown in light red, and the resulting elliptical fits to these 2D slices are in dark red. (d) Surface representations of the PAG (dark gray) enclosing S 2 (CA; light gray), shown with S 1 (brainstem surface; light gray. Three orthogonal slices of one subject's T1 anatomy volumes are shown. Light red depicts the actual PAG contour at p = .96 (S ). Dark red depicts the boundary of the estimated PAG segmentation (S ) reconstructed from transforming the 2D elliptical fits back to their 3D locations. (e) Left: Mean ± SD of elliptical parameters (a = left–right semi‐axis; b = dorsal‐ventral semi‐axis; ε = aspect ratio) across subjects. Right: Mean ± SD elliptical parameters a (orange) and b (green) along the length of the cerebral aqueduct, normalized to the mean length of CA

A second 3D representation of the PAG was also obtained as a smooth iso‐probability surface, S , of the gray‐matter probability map at using the methods described above (Figure 6d). Upon grossly comparing S and S , the former was, as expected, very similar to the elliptical fit described above.

Measurements of collicular thickness

We used path distances along the streamlines to create a new distance map estimating physical tissue depth, obtained by contour integration along each streamline. The path integrals, which oversample the volume, are then re‐gridded using a Delaunay triangulation approach (Amidror, 2002) onto the anatomical reference volume (Figure 2d). This level‐set depth stopped increasing at CA, whereas Euclidean depth (Figure 2a) continued to increase into the center of BS.

For the 12 subjects from which we collected gray‐matter probability maps to delineate the PAG, we measured collicular tissue thickness using three different depth metrics: Euclidean distance, level‐set path distance, and normalized w distance. Then, collicular thickness was calculated by integrating path distance along each streamline from the surface of the colliculi to S (7A). Table 1 displays mean collicular thickness values at six particular locations: the four peaks of left and right SC and IC, and the two sulci between the left and right colliculi of each SC/IC pair. Our method for automatically selecting these points is as follows (Figure 7b).

TABLE 1

Mean ± SD for tissue thickness values in left (L) and right (R) colliculi and intercollicular sulci (ICS) for SC and IC in all 12 subjects using three different depth metrics

	SC			IC
Units	L	R	ICS	L	R	ICS
Physical (Euclidean)	2.64 ± 1.59	3.09 ± 1.83	1.04 ± 0.70	3.35 ± 0.48	3.11 ± 1.17	0.88 ± 0.35
Physical (streamline)	4.19 ± 2.84	4.80 ± 2.79	1.23 ± 0.66	3.18 ± 0.32	2.75 ± 1.04	0.93 ± 0.43
Normalized (w)	0.43 ± 0.16	0.42 ± 0.13	0.26 ± 0.13	0.39 ± 0.04	0.32 ± 0.12	0.20 ± 0.07

FIGURE 7

(a) Surface representations of brainstem (S 1) of four subjects scanned at 9.4T, overlaid with collicular tissue thickness maps. (b) Determination of collicular peaks and intercollicular sulci. The collicular surface S 1A was divided into four quadrants using a quasi‐sagittal (blue) and axial (fuchsia) plane. A base plane (yellow) through the vertices with zero curvature was defined for each colliculus, and collicular peaks (red) were defined to be the vertices with greatest distance to each base plane Intercollicular sulci (cyan) were vertices on the mid‐sagittal plane closest to the line between left and right peaks

First, for each subject, we isolated the vertices of S 1 of corresponding to the colliculi. We manually created a third segmentation label for the voxels of SC and IC in ITK‐SNAP that grossly segregated the surfaces of all four colliculi. We then mapped these voxels by nearest‐neighbor association to the BS surface. For completeness, we then expanded this surface by an additional 2‐mm manifold distance, thus creating S 1.

Second, S 1 was used to locate the peaks of each colliculus by dividing its vertices into four quadrants, each containing one colliculus, based on a set of planes. A quasi‐sagittal plane separating left and right colliculi was specified as the plane of best fit through the vertices of the CA surface model, S 2. To calculate an axial plane separating SC and IC, the collicular surface vertices were divided into bins of 0.5‐mm width according to their axial position. The midpoint of the bin with the most negative mean curvature defined the axial plane of separation. We used these planes to divide S 1 into four subsets, each containing a single collicular surface. Next, we defined a base‐plane for each colliculus by taking all vertices with positive curvature values and fit these with a plane using a least‐squares method. We then defined the peak location as the point with maximum Euclidean distance from this base‐plane. We also experimented with defining the collicular peaks as the vertices with greatest positive curvature but found that these points of maximum curvature varied significantly across individual colliculi as well as across subjects. In contrast, the peak distance definition was far more stable, and was therefore used as our standard approach.

Third, for the intercollicular sulcus (ICS) of SC, we calculated the line between the peaks of left and right SC calculated as described above. The intersection point between this line with the plane‐of‐best‐fit through the CA was located, and the vertex of S 1 closest to this point, within 0.25 mm of the plane of best fit, was defined to be the center of the intercollicular sulcus. For the inferior sulcus, we repeated the same process with left and right peaks of IC.

Once we obtained the six points of interest, we measured each of their thickness values using three different metrics. For Euclidean distance, we used the distance along the surface normal from S 1 to S PAG. For level‐set physical distance, we used integrated path distance along the streamline originating from each point to the point where it intersected the PAG. For normalized distance, we linearly interpolated the w‐distance map at the point where the streamlines of interest intersected the PAG. Lastly, we applied a Gaussian‐weighted smoothing kernel (2‐mm diameter, 1‐mm full‐width‐at‐half‐maximum [FWHM]) to each of these points to obtain estimates of the six thickness values. For each of n vertices within 1 mm of the peak or sulcus in question, its thickness value t was given a weight A according to its manifold distance d to the center point of the kernel: , where is the standard deviation of the Gaussian function, and FWHM . The final thickness value was then calculated as the weighted sum, .

RESULTS

Euclidean versus level‐set methods

Unlike the kernels generated using the cylinder method, the level‐set streamlines avoided overlap through their nonlinearity and compression with increased depth (Figure 4a,b). The streamlines in the colliculi exhibited various degrees of non‐linear deformation depending on starting location and their depth into BS. With respect to superficial location, the greatest nonlinearities were expected to occur in streamlines that originate from regions where the curvature varies rapidly (Figure 4a,b). Low amounts of nonlinearity at the lateral edges of collicular tissue were observed, gradually increasing inwards towards the concave tissue between the colliculi.

With respect to metrics of physical tissue depth, the relationship between Euclidean and level‐set depth was sub‐linear (Figure 4d), reflecting increasing curvature in deeper collicular tissue. Additionally, level‐set depth, or path distance along the streamlines, stopped increasing at CA, whereas Euclidean depth continued to increase into the center of BS (Figure 2a,d).

We established a metric Δd to quantify the deviation between level‐set and cylinder depths (Figure 4c). In addition to measuring deviation between the two methods, this comparison also quantified the non‐linearity of the level‐set depth trajectories. At each point along each streamline, Δd was the orthogonal distance to the S 1 normal vector of the vertex from which the streamline originated. Greater values of Δd at increasing level‐set depth were consistent with the warping observed above. Δd was mapped back to S 1 at various depths (Figure 4e). Non‐linearity was very small (<<1 mm) for depths <3 mm and became substantial (>1 mm), at ≥5 mm depth. Largest error (and non‐linearity) was in intercollicular regions, demonstrating where the Euclidean cylinder method is least appropriate.

Depth profiles of visual stimulation‐ and saccade‐evoked activity in SC

In superficial SC, our level‐set depth was very similar to the original Euclidean depth. Our measurements of the BOLD response as a function of depth confirmed similar depth profiles, including the previous finding that the mean depth profiles across subjects for saccadic activity was shifted deeper into SC compared with those obtained from visual stimulation. This is consistent with the functional organization of SC delineated in animal models (Fuchs, Kaneko, & Scudder, 1985; Sparks & Hartwich‐Young, 1989; Wurtz & Albano, 1980). The profiles generated using our new method did not show significant differences from those utilizing the cylinder method (Figure 5). The laminar depth profiles generated using our streamline kernels and normalized w depth also demonstrate the depth differences between stimulation‐ and saccade‐evoked activity to a similar level of significance as the other methods.

Characterization of the PAG

Upon examination the gray matter probability maps of the PAG, we found that a large proportion of the dorsal midbrain is attributed a high probability value, even including superficial tissue that would be considered to be colliculus (Figure 6a). However, after using a combination of the streamline depth kernels and Euclidean distance from CA, we were able to localize the PAG within this ambiguous region with a relatively small degree of uncertainty (Figure 6b). The PAG for all subjects was generally able to be fitted with ellipses resulting in a root‐mean‐squared error of 0.65 mm (Figure 6c), which did not show any significant variations within the slices of individual subjects. Ellipse parameters a, b, and ε (Figure 6e) had coefficients‐of‐variation 0.37 ± 0.08, 0.15 ± 0.06, and 0.49 ± 0.34, respectively, demonstrating that the elliptical fits were significantly more stable across subjects in the dorsal‐ventral direction compared with the left–right direction. Rotation parameters were negligible. Visualization of the PAG contour and its associated surface, S , reveals that this nucleus is not necessarily a uniform tube surrounding CA. Rather, its thickness, particularly along the left–right axis, varies along the length of CA (Figure 6d,e).

The vertices of the PAG contour (S PAG) and our estimation of the contour boundary using ellipses (S ellipse) were very similar. Mean vertex separation, averaged across subjects, was 0.46 ± 0.18 mm. A few vertices were widely separated; maximum separations were 7.2 ± 5.0 mm. However, the fraction of these outlier vertices, with separations >2 mm was very small, mean across subjects was 4.3 ± 2.7%. Thus, the PAG is well characterized as a series of smooth ellipses that vary along the length of the CA. Our elliptical fitting procedure smoothed out many of the irregularities present in the actual PAG contours derived from the gray‐matter probability maps (Figure 6d), though some inhomogeneities are still observable in orthogonal cross‐sections of S ellipse. However, particularly visible in a sagittal view, these “dimples” occur in the plane of the 2D ellipses, or in a direction that is normal to CA, reflecting some natural variance in the elliptical parameters along the length of CA. Overall, our results confirmed the utility of Euclidean distance from CA as a reference depth metric for this region of the midbrain.

Measurements of collicular tissue thickness

Thickness maps calculated using level‐set physical distance from S 1 to S showed that the tissue of the colliculi tends to be thickest at its crown, as expected, but also extending laterally, particularly the superior lateral edges for SC. The tissue inbetween left and right colliculi, the intercollicular sulcus, was thinnest (Figure 7a). Mean data across subjects for tissue thickness at the four peaks of SC and IC as well as the two sulci between left and right SC and IC are shown in Table 1. Interestingly, the variability of SC thickness was 2.77 mm of level‐set depth, while IC thickness was significantly more stable at 0.97 mm. A similar pattern was observed for superior and inferior intercollicular sulci, which had variabilities of 0.88 and 0.46 mm, respectively. For absolute thickness measures, there were no significant left–right thickness asymmetries. Thickness values at the intercollicular sulci were very similar for inferior and superior colliculi. The normalized thickness metrics showed lower coefficients of variation in general (Table 2). Left–right symmetry was also more clearly suggested by normalized thickness.

TABLE 2

Coefficient of variation (COV) of tissue thickness values using three different depth metrics across all 12 subjects

	SC			IC
Units	L	R	ICS	L	R	ICS
Physical (Euclidean)	0.602	0.593	0.671	0.144	0.377	0.397
Physical (streamline)	0.676	0.582	0.539	0.101	0.381	0.465
Normalized (w)	0.364	0.298	0.513	0.094	0.369	0.342

DISCUSSION

We compared an algebraic level‐set scheme utilizing a normalized depth coordinate, w, to a Euclidean depth‐mapping scheme for midbrain tissue between the surfaces of the brainstem and the cerebral aqueduct in humans. Previously, we used computational cylinders and Euclidean nearest‐neighbor definitions of depth to estimate depth profiles of functional activity evoked by visual stimulation and attention (Katyal et al., 2010; Katyal & Ress, 2014; Savjani et al., 2018). However, the cylinder method was limited by contention of depth kernels, rendering it susceptible to oversampling deep tissue within individual colliculi and undersampling tissue in the inter‐collicular sulcus. In contrast, the resulting nonlinear coordinate system in native brain space using our level‐set method is specifically adapted to the 3D shape of the tectum. Our novel streamline depth kernels never crossed paths by nature of the pseudopotential and compress or expand along their depth, which avoided mis‐association between deep and superficial voxels using the cylinder method. The streamlines give a one‐to‐one association between the brainstem surface and the cerebral aqueduct (CA) at all spatial locations.

We also defined a new metric for physical tissue depth using streamline path distance. Both this level‐set physical depth and the normalized w depth provide an analytic depth coordinate for each voxel. Additionally, the streamline depth kernels enable surface‐based smoothing of the functional activity, provide mappings between the superficial surface (S 1) and deeper tissue, and allow quantitative data, such as gray matter probability, at different depths to be projected onto the surface (Dale, Fischl, & Sereno, 1999; Glasser et al., 2016). The parcellation of tissue into kernels provided depth relationships between deep and superficial tissue at a single spatial location, or over any functionally defined superficial surface. For example, the retinotopic organization of superficial superior colliculus (SC) (Cynader & Berman, 1972) could be related to the putative somatotopic organization in the deep layers of human SC suggested by experiments on animal models (Clemo & Stein, 1991; Jay & Sparks, 1987; Nagy, Kruse, Rottmann, Dannenberg, & Hoffmann, 2006; Wallace, Wilkinson, & Stein, 1996).

Overall, the streamline kernels are superior to the cylindrical kernels in their ability to create logical associations between tissue. However, we show that both the Euclidean and level‐set depth metrics have their merits in appropriate regions of the midbrain. Usage of streamline depth kernels and level‐set physical depth produced depth profiles of visual stimulation‐ and saccade‐evoked activity in superficial SC similar to those previous obtained using cylindrical kernels and Euclidean distance from the surface of the brainstem, d 1, as the reference depth metric (Figure 5). This similarity was maintained throughout the depth range over which we calculated the depth profiles, 0–3 mm, and negligible differences between the methods in superficial tissue suggest that either method works well at these superficial depths. Meaningful functional activity did not extend deep enough into SC for us to determine any significant differences between Euclidean and level‐set depth metrics beyond these depths. This is consistent with our findings in Figure 4d, which demonstrate that substantial differences (e.g., >0.5 mm) between Euclidean and level‐set path distance metrics only emerge at depths >3 mm. Depth profiles generated using normalized w distance and streamline kernels also produced results with similar degrees of significance and variability compared to those generated with physical measures of depth. According to Figure 4d, the medial regions of the colliculi are most susceptible to significant deviations between cylindrical and streamline sampling methods. Notably, upon visual inspection, the elliptical ROIs obtained in our fMRI experiments typically did not cover medial SC, which may explain the lack of significant differences in the depth profiles obtained using the two methods. On the other hand, the anatomical structure of the PAG is best defined using depth values according to Euclidean distance from CA, or d 2 (Figure 6). As d 2 is straightforwardly defined using the boundary of CA as a reference surface, this is in agreement with the qualitative description of the PAG as roughly a cylindrical tube surrounding CA. Thus, Euclidean distances from an appropriate reference surface are sufficient to delineate both superficial colliculus and the deep PAG.

Because of the differences between the methods that emerge at greater depths within the midbrain, characterized by our metric Δd, we predict that our streamlines and level‐set depth metrics will provide the greatest utility in regions of intermediate depth, namely deep SC and IC. This hypothesis is made on the basis that our level‐set streamline kernels more accurately follow the natural curvature of midbrain tissue. Particularly, Figure 4d demonstrates that the greatest deviations between Euclidean and level‐set methods are located in the medial regions of the colliculi, where the surface normal of the brainstem tissue point away from CA, indicating where use of cylindrical depth kernels may be least appropriate. In regard to the laminar structure of the colliculi, isosurfaces of our normalized w distance may provide a more accurate representation of the anatomical morphology of midbrain tissue compared with Euclidean definitions of tissue depth. This should be validated through comparison studies between high‐resolution in vivo brain imaging and histology on post‐mortem brains, similar to (Loureiro et al., 2018). Further confirmation of the utility of our level‐set method in deeper midbrain structures, particularly for fMRI experiments, should be achieved through studies such as those involving reach‐related activation in SC (Linzenbold & Himmelbach, 2012; Liu, Duggan, Salt, & Cordeiro, 2011; Sparks & Hartwich‐Young, 1989). Given the ability of our w‐depth profiles to differentiate laminar functional activity in superficial SC, we believe that the applications of our method in deeper tissue are promising.

We measured collicular thickness using our level‐set definition of depth. While thickness measured in physical units (mm) varied between individuals, measurements in units of normalized depth demonstrated the least variability, testifying to the normalization that this metric provides between brains of different sizes. SC tended to have similar thickness, but greater variability compared to IC. While the mean thicknesses of left and right colliculi were not significantly different, left–right symmetry varied among individual subjects. Some subjects had differences in thickness greater than 2 mm, while others had differences <0.25 mm. The left–right differences in thickness were somewhat smaller for SC than IC.

However, from visual inspection of the PAG contour maps within the context of each subject's anatomical volume, left–right asymmetry was not apparent. This disparity between expectation and qualitative observation, and our data are subject to two interpretations. First, what appears to be qualitatively symmetrical might not be in terms of path distance along nonlinear streamlines, which is not immediately evident purely from visual inspection. Second, and more importantly, the asymmetry tended to be an effect of which point on the BS surface was chosen to be the “peak” of each colliculus. Our method of choosing the point with the furthest straight‐line distance to the base plane of each colliculus has an inherently Euclidean basis. Alternatively, taking the point with the greatest streamline distance to the base plane would have introduced a bias towards the very method we are attempting to test. It would also favor choosing points near the lateral edges where the streamlines are least perpendicular to the base plane. The desired output of the peak‐finding algorithm is unclear at this point, especially when we have yet to fully characterize the location and lateral boundaries of the colliculi.

To the best of our knowledge, this work is the first to assess human collicular thicknesses in vivo, and these measurements, expanded to a larger pool of healthy brains, may be useful for gross delineation of neurodegenerative pathology that affect the colliculi such as glaucoma or aging (Crish, Sappington, Inman, Horner, & Calkins, 2010; Liu et al., 2011; Ouda & Syka, 2012). Though our data may appear to imply that some of our subjects suffered from unilateral degeneration, our method is one of many possibilities. Developing a way to accurately measure collicular tissue thickness depends on our ability to localize the entire boundary around these nuclei. Again, comparison studies between brain imaging and histology on post‐mortem brains (Loureiro et al., 2018) will be necessary to provide a complete analysis of collicular topology.

We strove for a method to relate observed function back to the anatomical structures from which they originate. Studying brain function within a structural context, as in a standardized atlas space (Evans et al., 1994; Mazziotta et al., 2001; Talairach & Tournoux, 1988), allows for localization of relevant structures and, consequently, generalizability across populations and studies. For example, the MNI atlas allows non‐linear warping of 3D MRI brain scans from individuals onto a common template and has become a standard tool for reporting fMRI results in volume space. The PAG particularly lends itself for atlas comparisons, since this brain structure is relatively stable across subjects and with age (Keuken et al., 2017). Our experiments assumed that the tissue probability maps we used were well‐suited to study the PAG. Usage of these maps has been validated using other deep brain nuclei including the thalamus, caudate, putamen, substantia nigra, subthalamic nucleus, red nucleus, and cerebellar dentate (Lorio et al., 2016). To our knowledge, this work is the first application of these methods to the PAG. Though these maps were not explicitly designed for localization of the PAG, small spatial uncertainty values (Figure 6b) suggest that this application is appropriate. Our more rigorous quantitative method confirms the use of high gray‐matter probability values for tissue classification. Furthermore, under the assumption that the PAG boundary can be fit by a continuous surface, our estimation scheme using ellipses worked well. We observed low root‐mean‐squared errors for individual elliptical fits, and low coefficients‐of‐variation for the ellipse parameters across subjects.

Further extension of our techniques to provide a fully 3D coordinate system may utilize normalized depth and angle from a point of origin on the mid‐sagittal plane through CA. This could potentially allow localization of midbrain nuclei such as the PAG using polar coordinate shapes such as dimpled or convex limaçons, which may provide a greater degree of accuracy than the simple ellipses we utilized in our work. Our efforts towards characterizing and subsequently reconstructing the PAG provide a first step in the direction towards creating a generalized atlas of midbrain nuclei (Figure 6). Using CA as a central axis, as we have done, could form the basis of a quantitative atlas of the midbrain in a normalized coordinate space such as the MNI‐305 (Mandal, Mahajan, & Dinov, 2012), but at higher spatial resolution. While we concentrated our evaluation of this depth‐mapping scheme to the dorsal portion of midbrain, our method can be also effective for ventral midbrain, for example, the red nucleus. It may be possible to interpolate the neuraxis between the exit of CA into the fourth ventricle and the geometric axis of the spinal cord at the base of the brainstem to facilitate similar methods in pons and medulla.

Our depth analysis assumed that the laminar structure of SC is well‐delineated by a depth metric derived from a combination of SC's superficial surface and the CA. In order to minimize oversampling and overlapping depth projections, selecting an appropriate reference structure was key for defining depth in the midbrain. The CA runs rostro‐caudally through the midbrain ventral to the colliculi. It emerges developmentally from the neural canal, the cavity in the center of the neural tube. The midbrain develops around the CA, thus positioning it as an important anatomical feature in the midbrain. Ideally, the resulting isosurfaces of w should be parallel to the laminae of the colliculi. However, our choice of surfaces to calculate depth coordinates was also a matter of convenience. Although our present results are consistent across individuals and demonstrate the utility of our depth mapping approach in vivo, it remains to be seen how well the depth‐mapping scheme aligns with laminarly organized collicular cytoarchitecture. In future work we, will apply our depth mapping approach to publicly available datasets, such as BigBrain and CIT168 MNI, and high‐resolution post mortem MRI and histology (Loureiro et al., 2018; Sitek et al., 2019), where individual collicular laminae are distinguishable. This will allow us to validate our method's sensitivity to ground truth depth‐dependent cytoarchitecture.

CONCLUSIONS

We established a series of morphologically reasonable methods for determining tissue depth in human SC and IC based solely on structural MRI. The level‐set scheme provides coordinates for both normalized and physical depth that accommodate the complex anatomical structure of midbrain tissue, and the associated streamline method was able to distinguish function in the superficial layers of SC at the same precision as our previous Euclidean method. Euclidean distance from CA provides an accurate way to localize the PAG. We also believe that our work is the first analytic method for in vivo determination of anatomical tissue thickness in human midbrain. Because the streamline depth kernels logically follow the curvature of the tissue to establish unique relationships between all layers, our method should be useful to describe associations between tissue at different depths, such as between the deep multi‐sensory layers of superior colliculus and the superficial retinotopically organized layers. Furthermore, depth‐analysis studies using our normalized depth coordinate applied to a large population of subjects may enable a fully 3D atlas of human brainstem.