A surface-in gradient of thalamic damage evolves in pediatric multiple sclerosis.

OBJECTIVE: Central nervous system pathology in multiple sclerosis includes both focal inflammatory perivascular injury and injury to superficial structures, including the subpial region of the cortex, which reportedly exhibits a gradient of damage from the surface inward. We assessed how early in the multiple sclerosis course a "surface-in" process of injury suggesting progressive biology may begin.
METHODS: We focused on the thalamus, which notably has both a cerebrospinal fluid (CSF) interface and a white matter interface. Thalamic volume trajectories were assessed in a prospectively followed cohort of children from initial presentation with either multiple sclerosis or monophasic acquired demyelination, and healthy controls. Voxelwise volume changes were calculated using deformation-based morphometry, and analyzed in relation to distance from the CSF interface by mixed effects modeling and semiparametric smoothing methods.
RESULTS: Twenty-seven children with multiple sclerosis and 73 children with monophasic demyelination were prospectively followed with yearly brain scans (mean follow-up = 4.6 years, standard deviation = 1.9). A total of 282 healthy children with serial scans were included as controls. Relative to healthy controls, children with multiple sclerosis and children with monophasic demyelination demonstrated volume loss in thalamic regions adjacent to the white matter. However, only children with multiple sclerosis exhibited an additional surface-in gradient of thalamic injury on the ventricular side, which was already notable in the first year of clinical disease (asymptote estimate = 3.01, 95% confidence interval [CI] = 1.44-4.58, p = 0.0002) and worsened over time (asymptote:time estimate = 0.33, 95% CI = 0.12-0.54, p = 0.0021).
INTERPRETATION: Our results suggest that a multiple sclerosis disease-specific surface-in process of damage can manifest at the earliest stages of the disease. ANN NEUROL 2019;85:340-351.

Multiple sclerosis (MS) is a chronic inflammatory disease of the central nervous system (CNS) and a major cause of neurological disability in young adults.1 Classical pathological hallmarks include multifocal perivascular demyelinating white matter (WM) lesions,2 which have been associated with relapsing–remitting disease biology. However, a more gradual and diffuse progressive injury process (referred to here as a "non‐relapsing progressive injury") is thought to be a major contributor to the accrual of irreversible neurological deficits.1, 3 Cortical and deep gray matter involvement is increasingly recognized as an important feature of this progressive injury.4, 5 Pathology studies of the subpial cortical injury describe demyelination, neuronal loss, microglial activation, and a relative paucity of inflammatory cells, as well as a pattern of damage that follows a "surface‐in" gradient.4, 5 These observations are supported by recent magnetization transfer ratio (MTR) imaging studies in MS reporting a graded surface‐in abnormality starting at the brain–cerebrospinal fluid (CSF) interfaces of both the cortical surface and the ventricular surfaces.6, 7, 8 Together with the correlation between the severity of subpial cortical damage and the degree of meningeal inflammation,9, 10, 11 these findings have raised the hypothesis that soluble factors, potentially released by immune cells in the meninges, diffuse into the CSF and into the superficial brain layers such as the subpial cortical tissue, contributing to the observed gradient of pathological changes.12, 13, 14, 15, 16

Regardless of the underlying mechanism, it is unknown when in the disease course such a nonrelapsing progressive biology of MS begins.17, 18, 19, 20, 21, 22, 23 Whereas clinical emergence of nonrelapsing disease progression tends to be appreciated later in the disease, the biological process or processes contributing to it are now thought to start earlier, supported by several studies of adult onset MS showing that reduced brain volumes and progressive brain atrophy are appreciated from initial clinical presentation.19, 20, 21, 24 Prior work in pediatric onset MS patients further supports this concept, with the demonstration of reduced brain volumes, diffusion tensor imaging evidence of disruption of the integrity of normal appearing tissue, and failure of normal brain growth, all observed starting at clinical disease onset in these very young patients.17, 18, 23

We hypothesized that a surface‐in process of damage adjacent to the CSF (as a possible reflection of the nonrelapsing progressive MS biology) occurs very early in the disease course and thus will already manifest in children presenting with their first clinical episode of MS. Because pathology work has suggested that the subpial cortical injury is a particular feature of MS (not found in a range of other inflammatory CNS conditions),25, 26 we further hypothesized that CSF surface‐in damage would be specific to MS and would not manifest in children with monophasic acquired demyelinating syndromes (monoADS) studied in parallel. We focused on the thalamus for several reasons: (1) it is one of the most atrophic structures in early adult onset MS (reviewed in Minagar et al27) and is the brain region with the greatest deviation from age‐expected growth in pediatric onset MS, with volumes that are abnormally low already at clinical onset, and with progressive atrophy that continues thereafter;17, 18, 28 (2) it has a CSF–interface on the ventricular side and a WM interface more laterally, providing a unique opportunity to examine whether distinct injury processes manifest at the two interfaces; and (3) its regular shape and reasonably smooth and well‐defined boundaries enable more precise volumetric analyses, including less confounding by partial volume effects (the mixing of ≥2 tissues, such as gray matter and CSF, in voxels along an interface), as compared to the cortex. To date, pathologic studies of the thalamus in MS have reported a relatively prominent involvement of the periventricular regions, significant neuronal loss, and a degree of inflammation intermediate between the inflammation observed in typical perivascular WM lesions and subpial cortical lesions. There have been no tissue‐based studies directly assessing whether the thalamus exhibits a gradient of injury. As noted above, MTR imaging studies in adult onset MS have reported a graded surface‐in abnormality both at the cortical surface and at the ventricular surface, which would be expected to include thalamic structures.6, 7, 8

Subjects and Methods

Participants

We studied 100 children enrolled as part of the Canadian Pediatric Demyelinating Disease Study (CPDDS; Table 1). All participants for this study were recruited at a single site (Hospital for Sick Children, Toronto, Ontario, Canada) and were required to have serial MRI studies that passed a predefined quality control process. All were recruited at the time of an incident clinical episode suggestive of ADS and then followed prospectively with clinical assessments and annual coregistered brain magnetic resonance imaging (MRI) scans to ascertain the diagnosis of MS (n = 27) versus monoADS (n = 73). MRI scans of 258 healthy control (HC) children from the National Institutes of Health Study of Normal Brain Development dataset (NIHPD HC) were included as a normative control group29 together with 24 healthy children with a negative history for neurologic, medical, or psychiatric diseases (CPDDS HC), imaged on the same scanner as the CPDDS patients at 2 time points, 12 to 24 months apart, as per our prior work confirming the absence of scanner‐related differences between the CPDDS and NIHPD cohorts.18 Guardians and participants provided written informed consent. Younger children provided verbal assent. The study was approved by ethics boards of participating institutions.

Table 1

Demographic and Clinical Characteristics of the Pediatric Study Cohort

Characteristic	MS	monoADS	HCs
Participants, n	27	73	258 NIHPD HCs, 24 CPDDS HCs
F, n (%)	18 (67%)	33 (45%)	160 (57%)
Age at baseline, yr, mean (SD), range	12.56 (2.83), 5.90–15.90	10.07 (3.60), 2.28–15.94	11.51 (4.00), 4.23–21.91
Age at fist scan, yr, mean (SD), range	12.63 (2.83), 5.92–15.91	10.32 (3.60), 3.84–16.17	11.51 (4.00), 4.24–21.91
Years of follow‐up, mean (SD), range	4.55 (2.04), 2.00–8.92	4.85 (1.88), 1.93–8.97	3.02 (1.08), 1.01–5.4
Scans, n	119	332	684
Number of scans per subject, mean (SD)	4.41 (1.93)	4.55 (1.68)	2.42 (0.50)
Time elapsed between scans, yr, mean (SD)	1.14 (0.60)	1.05 (0.44)	1.98 (0.87)
Baseline whole thalamic volume normalized to brain, mm3 (SD)	17,384.4 (964.4)	17,290.3 (1,055.3)	17,490.8 (1,088.8)
Whole thalamic volume change, % per yr	−0.971	0.000	0.073

CPDDS = Canadian Pediatric Demyelinating Disease Study; F = female; HC = healthy control; monoADS = monophasic acquired demyelinating syndrome; MS = multiple sclerosis; NIHPD = National Institutes of Health Study of Normal Brain Development; SD = standard deviation.

MRI Protocol

All MS, monoADS, and HC participants in the CPDDS study were scanned on a 1.5 T Twinspeed Excite 12.0 scanner (GE Healthcare, Waukesha, WI). The protocol included a 3‐dimensional (3D) T1‐weighted (T1w) radio frequency (RF)‐spoiled gradient‐recalled echo (repetition time [TR] = 22 milliseconds, echo time [TE] = 8 milliseconds, 30° flip angle, 0.98 × 0.98 × 1.5 voxel size). The NIHPD HCs were scanned at 6 centers on GE and Siemens (Erlangen, Germany) 1.5 T scanners. All acquisitions included a 3D T1w RF‐spoiled gradient‐recalled echo (TR = 22–25 milliseconds, TE = 10–11 milliseconds, 30° flip angle, 1 × 1 × 1 voxel size).

Image Processing

Thalami were segmented on the T1w images of each participant at baseline and at every yearly follow‐up, with atlas‐based segmentation using Automatic Nonlinear Image Matching and Anatomical Labeling (ANIMAL) from the MINC toolkit,30 which would serve to limit possible bias that may be introduced due to the differences in boundary contrast. Images were cropped to a region of interest that included both thalami and the surrounding tissue in a rectangular prism extending 10 voxels on all 6 sides from the maximum extent of the thalamus. Using fine nonlinear registration (ANTs toolkit, SyN algorithm, 200 × 200 × 100 × 50 iterations, cross‐correlation, 1 mm smoothing of deformation field, 2 mm smoothing of gradients; registration generally ran to convergence at each level), the cropped T1w images were nonlinearly registered to the baseline image and deformation‐based morphometry (DBM) was applied to calculate the volume change of each voxel after the registration (the local Jacobian determinant). An advantage of DBM31 is that the thalamus is only segmented once for each subject, such that changes over time are more likely to reflect actual changes in volume rather than the potential influence of changes in contrast of any given thalamic border. Registered images were averaged for each participant, producing subject‐specific templates, which were then registered nonlinearly to the Montreal Neurological Institute (MNI) pediatric template and averaged to define a grand average template. Each transformed image was inspected to ensure that the registration was of high quality. Individual Jacobian determinant maps were also transformed into grand average space, along with thalamus and CSF masks. This procedure gives the option of comparing groups either in the common grand average space, as is conventionally done, or in their individual subject‐specific spaces. Whereas the former procedure aligns voxels with similar positions relative to surrounding anatomy, analysis in the subject‐specific space has the added value that it preserves true distances, measured in millimeters. The comparison of the group trajectories in the common and the subject‐specific space in our study revealed no substantial differences. Because we hypothesized a process mediated by factors diffusing from the CSF, we considered the absolute distance from the surface as a more relevant factor, and present our analyses on the data obtained in the subject‐specific space.

Analytic Approach

We first compared whole‐thalamic volume change between children with MS and those with monophasic demyelination relative to age‐ and sex‐matched HCs. Whole‐thalamic volume measurements were made on each scan using the thalamic masks, and normalized to whole‐brain volume, as per our previous work.18 To account for differences in age and sex that could influence the growth trajectory, we performed nonparametric normalization of the two patient cohorts by subtracting from each participant's measurements the average of all available HCs of the same sex and age ± 6 months at the time of scanning. The resulting normalized values were analyzed using a linear mixed effect model, with fixed effects for patient group (MS or monoADS), elapsed time (the time between first clinical presentation and the present scan), and interactions between these, and a random intercept for each subject.

Data generated using the DBM technique31 provide measurements of local quantitative tissue expansion or contraction for each voxel. Most DBM analyses employ comparisons at the voxel level, which produces statistical maps, but requires a large number of statistical comparisons and thus the use of a method for false‐positive control. Instead, we found the average deformation in the thalamus at fixed distances from CSF, then analyzed these data with diffusion‐inspired mixed effects models. This approach allowed us to rigorously test our hypothesis with single models, limiting multiple comparisons and obviating the need for qualitative interpretation of statistical maps.

A Euclidean distance transform was applied to the CSF masks, producing a voxelwise map of distance from the ventricular surface, from which the thalamus was divided into 1 mm‐thick concentric surfaces of equal distance from CSF (Fig 1). Because our registration employed a 1 mm smoothing kernel, and the Jacobian in a voxel is itself affected by its neighbors, measurements within 2 mm of CSF were excluded from analysis to avoid potential contamination by partial volume effects. Voxels that were > 12 mm from CSF were also excluded, as few participants had thalami large enough to provide data past this point.

Figure 1

Measurement of the distance from cerebrospinal fluid (CSF). Distance of each thalamic voxel is shown in millimeters from CSF, obtained after applying the Euclidean distance transform to the CSF mask.

The voxelwise volume changes were then averaged over each isodistant surface to obtain the average change in volume for tissue at each millimeter, moving in depth away from the ventricular surface. As was performed for the whole thalamic volume calculations, we accounted for differences in age and sex that could influence the growth trajectory, using nonparametric normalization of the patient measurements to HCs by subtracting from each participant's measurements the average of HCs of the same sex and age at scan ±6 months. Through this approach, we also accounted for the intrinsic anatomical heterogeneity of the thalamus and isolated the abnormal volume changes in different areas of the thalamus from expected regional differences in normal growth.

We then considered whether thalamic growth in our patient groups differed from normal thalamic growth by first evaluating thalamic tissue close to, and furthest from, CSF. Specifically, volume changes for each participant at each time point were averaged in thalamic tissue 2 to 6 mm from CSF (CSF side) compared to tissue 8 to 12 mm distant from CSF (WM side). These simple summary measurements were compared using a linear mixed effect model with fixed terms for side (CSF or WM), patient group (MS or monoADS), and their interaction, plus a random intercept for each participant.

To isolate the MS‐specific pattern of abnormal growth, we again performed nonparametric normalization of the local volume change data, this time subtracting from the MS participants the average values of monoADS patients of the same sex and age at scan ±6 months (matching as many monoADS patients as available for each MS participant).

We considered that the graded (surface‐in) MS‐specific pattern of the injury observed at the CSF surface of the thalamus might be mediated by diffusion of one or more toxic factors from the CSF, with higher concentrations of such factors associated with more tissue damage. Such a process would be expected to follow Fick's second law of diffusion.32 To model this, we utilized a solution to Fick's second law, describing the concentration of a diffusible substance away from a surface of constant concentration, as illustrated in Figure 2. When the concentration at the source of a diffusible factor remains constant (as may be the case with CSF acting as a replenishable reservoir), the solution to Fick's second law takes the form of the error function, erf, which has the general shape of a sigmoid.

Figure 2

Example of Fick's law of diffusion. Example shows the form of Fick's law of diffusion. The concentration of the diffusible molecule is initially high and drops with distance from a large reservoir, creating a sigmoidal curve (orange). With a negative coupling between concentration and toxic effect (eg, volume change), the local effect mirrors the concentration (blue). The parameters describing a sigmoid function are: a, the asymptote, corresponding to the difference between the ceiling and the bottom value; b, the steepness or slope of the sigmoid; c, the crossover point where the sigmoid's slope begins decreasing. Curves plotted for a = 4, b = 0.35, and c = 6.5 are shown.

Our model was:where ΔV is the change in volume, I is an intercept term corresponding to global thalamic volume change not associated with distance from CSF, a is the lower asymptote, corresponding to the distance‐related volume change at the ventricular surface, b is the steepness or slope of the sigmoid, and c is the crossover point where the sigmoid's slope begins decreasing and distance is the distance from the nearest CSF, in millimeters. Unknown parameters, such as the concentration of the putative soluble toxic factor and the concentration‐related toxicity to tissue, are not included in this model.

We used Equation 1 in nonlinear mixed effects (NLME) models fit to the MS‐specific local volume changes. The statistical model:where addition of γ provides the subject‐specific random intercept and ε is an error term, tests how the distance from CSF affects the change in volume in participants who were diagnosed with MS, relative to those with monoADS. To investigate how the change in volume varied with time after disease presentation, we generated a second model including the interaction of time with the asymptote a of the sigmoid. We detected a nonlinear rate of atrophy, and therefore included a quadratic time–asymptote interaction, obtaining the following formula:where timeC is the time elapsed from disease presentation, centered by subtracting 4.55 years (corresponding to the mean duration of follow‐up). To compute the predicted change in volume over the first year, we calculated a contrast with timeC = −3.55 years. This analysis was then repeated using the date of MS diagnosis rather than the date of the first MS attack.

To further investigate the fit of the diffusion model, we devised an analysis designed to assess the pattern of any residuals (ie, data that did not align with a diffusible factor). To do so, we fit a generalized additive mixed effects model (GAMM) that assessed whether there was a systematic association between distance from CSF and the residual values from our diffusion model. Such an association would indicate the presence of effects that cannot be explained by diffusion. The GAMM had a random intercept for each participant and a fixed‐effect thin‐plate regression spline relating distance to residual value. We used a parametric bootstrap to compare this model against a null model with only a random intercept, which corresponds to the hypothesis that the variation is completely explained by the diffusion model. Finally, to ensure that the diffusion model was replicable using a different modeling technique, we fit a nonparametric model (LOESS) relating distance from CSF to local thalamic volume change.

We considered that our volumetric measurements could be artifactually affected by differences in the sharpness of the interfaces between the thalamus and the CSF, versus the thalamus and the WM. To address this in the analysis, we removed potential interface‐specific differences by replacing all voxels outside the thalamus in the cropped image of each subject with the intensity of mean normal‐appearing WM for that subject, then repeated the complete procedure, from nonlinear registration through statistical analysis, on these modified data.

The proportion of variance explained by the models, including random effects, was evaluated using the conditional R.2, 33 Data processing and statistical analyses were computed using Python (v2.7.11) and the LME4,34 NLME,35 GAMM,36 and PBKRTest37 packages in R (v3.1.3).

Results

A Pediatric MS‐Specific Gradient of Thalamic Injury as a Function of Distance from the Ventricular Surface

The basic demographic features of our study cohort are summarized in Table 1. Our study included 27 subjects with MS, 73 with monoADS, and a total of 282 HCs. We confirmed our prior findings of a negative impact on thalamic growth in both pediatric MS (p < 0.0001) and monophasic acquired demyelination (p < 0.0001) relative to HCs (data not shown). Figure 3A shows the relationship between distance from CSF and change in thalamic volume in each of the 3 cohorts of participants (MS, monoADS, HCs). Only in HCs did thalamic regions exhibit the expected tissue growth (positive values), most evident on the WM side of the thalamus. To define abnormalities in growth trajectories of the MS and monoADS cohorts, we normalized their results based on growth trajectories of age‐ and sex‐matched HCs (see Fig 3B). In these normalized data (which account for expected age‐ and sex‐related growth of different thalamic regions), values of zero indicate normal growth trajectories, whereas values different from zero indicate abnormal volume changes in MS or monoADS children. We found that children with monoADS exhibited failure of expected thalamic growth in the regions adjacent to the WM, whereas on the ventricular side there was no suggestion of abnormality. In contrast to children with monoADS, whose thalami were abnormal only on the WM side, children with MS exhibited 2 distinct thalamic abnormalities. The first was an abnormality on the WM side of the thalamus that was indistinguishable from the abnormality noted in the monoADS group. The second abnormality, unique to children with MS, was a gradient of atrophy that was most severe subjacent to the ventricular (CSF) surface. Specifically, as summarized in Figure 3C, decreased thalamic growth (relative to normal expected growth) was seen on the WM side of the thalamus in both children with MS (−2.9%, 95% confidence interval [CI] = −4.6 to −0.39, p = 0.0008) and children with monoADS (−2.3%, 95% CI = −3.3 to −1.4, p < 0.0001) and was of a similar magnitude (monoADS vs MS, p = 0.87). The MS group also demonstrated marked failure of normal growth of thalamic tissue closest to the CSF (MS: −4.7%, 95% CI = −6.7 to −2.5, p < 0.0001), which was not observed in the monoADS patients (monoADS: −0.48%, 95% CI = −1.5 to 0.65, p = 0.37; monoADS vs MS, p < 0.0001).

Figure 3

Changes in thalamic volume (%) compared to baseline as a function of distance from cerebrospinal fluid (CSF). (A) Trajectories for multiple sclerosis (MS; blue), monophasic acquired demyelinating syndromes (monoADS; green), and healthy control (yellow) cohorts. Shaded areas are the 95% confidence intervals (CIs) for the point estimates. Negative values indicate tissue loss, and positive values indicate tissue growth. Only the healthy control cohort exhibits thalamic growth. (B) Abnormal trajectories for MS and monoADS cohorts after normalization to healthy controls, where the horizontal line passing through 0 reflects the expected growth trajectory for healthy individuals. Both MS and monoADS cohorts fail to follow the expected growth in the more lateral thalamic regions adjacent to the white matter (WM), whereas only the MS cohort exhibits significant atrophy in proximity to the ventricular (CSF) surface. (C) Linear model results of the comparison between the average volume changes on the CSF side (2–6 mm from ventricular surface) and WM side (8–12 mm from ventricular surface) between the MS and monoADS groups. Both the MS and monoADS groups showed significant difference from normal growth on the WM side, whereas only the MS group showed a significant loss of volume on the CSF side. Error bars represent 95% CI.

Modeling the MS‐Specific Gradient of Thalamic Injury Supports a Diffusion‐Based Mechanism of Injury

We next subtracted the monoADS volume changes from volume changes of age‐ and sex‐matched children with MS, to isolate the MS‐specific pattern of thalamic abnormality (Fig 4A; green line). This highlighted an MS‐specific abnormality on the ventricular side of the thalamus, which was greatest in proximity to the ventricular surface and diminished in a graded fashion with increasing distance from that surface. We noted that a sigmoid NLME model fit the MS‐specific injury raw data remarkably well (see Fig 4A, blue line), which is consistent with a diffusion‐based mechanism of injury emanating from the CSF interface. Of the 4 parameters in Table 2 that describe a sigmoid function (intercept, asymptote, slope, and crossover point), the asymptote captures the distance between the floor and the ceiling values, therefore providing the amount of MS‐specific abnormal change that is described by the diffusion model. Though R values must be interpreted cautiously in nonlinear models, the sigmoid function R 2 was 0.71, indicating a good fit to the data. The GAMM model confirmed that the diffusion model fit remarkably well, identifying no apparent residual association between volumetric changes and CSF boundary distance, after accounting for the hypothesized diffusion‐related process (p > 0.9; see Fig 4B). The nonparametric LOESS fit revealed a similar sigmoidal pattern of greater volume loss in proximity to CSF unique to the children with MS (data not shown).

Figure 4

Fitted model for the multiple sclerosis (MS)‐specific pattern of damage. (A) The MS‐specific pattern of thalamic abnormality (green line) was calculated as the difference between the normalized volume trajectories of MS and monophasic acquired demyelinating syndromes (monoADS), and is shown as a function of distance from cerebrospinal fluid (CSF), with the shaded area representing the point estimate 95% confidence interval. Superimposed is the mixed effects model fit (blue line). The two lines closely overlap. (B) Results of the generalized additive mixed effects model fit (black) to the residuals (red) from the mixed effect model in A, using penalized splines. There is no evidence of residual association between volumetric changes and distance from CSF–interface after accounting for the hypothesized diffusion‐related process. (C) The mixed model estimates for the MS‐specific change in volume at 1, 2, 3, and 4 years from presentation (blue, green, red and purple lines, respectively). The MS‐specific gradient of thalamic atrophy is already detectable after 1 year of follow‐up and increases over time.

Table 2

Parameter Estimates of the Nonlinear Mixed Effects Model Analyzing the Change in Thalamic Volume in MS Patients Compared to Children with Monophasic Demyelinating Disorders, as a Function of the Distance from the Ventricular Surface

Parameter	Value	Standard Error	df	t	p
Intercept, I	3.019	1.91	982	1.58	0.11
Asymptote, a	4.10	0.10	982	4.12	<0.0001
Slope, b	0.29	0.12	982	2.41	0.016
Crossover point, c	8.48	0.79	982	10.68	<0.0001

Random effect: subject standard deviation = 4.82; residual standard deviation = 7.43; R 2 = 0.709. Values are compared to zero, which indicates no difference between MS and monophasic disease cohorts, both normalized for matched healthy controls. The intercept (I) represents the overall uniform atrophy; the asymptote (a) is the difference between the bottom and the ceiling values of the sigmoid, providing the amount of MS‐specific abnormal change in volume observed in proximity to the CSF interface; the slope (b) and crossover point (c) describe the shape of the sigmoid, and are related to the rate of diffusion in a diffusion model. The significantly positive value for the asymptote indicates significantly more atrophy in proximity to the CSF surface in MS compared to monophasic children.

CSF = cerebrospinal fluid; MS = multiple sclerosis.

Effect of Time Elapsed from Clinical Onset on MS‐Specific Changes in Thalamic Volume

Figure 4C and Table 3 show the results of NLME modeling including "time‐elapsed" since initial clinical presentation as an additional factor to explain the changes in volume (conditional R 2 = 0.71). The distance‐related change in thalamic volume subjacent to the CSF surface increased with time in the children with MS (0.33, 95% CI = 0.12–0.54, p = 0.0021; quadratic effect = 0.01, 95% CI = −0.088 to 0.11, p = 0.88; see Table 3), with evolving tissue loss in proximity to the ventricular surface already detectable in the first year of clinical disease (3.01, 95% CI = 1.44–4.58, p = 0.0002; Table 4). The alternative NLME model including the effect of time since diagnosis revealed only minimal differences compared to the one aligning patients at the time of clinical presentation (Supplementary Table). This was largely because a great proportion of MS patients (19/27, 70%) already met the 2017 McDonald diagnostic criteria for definite MS at the time of presentation, whereas the remaining patients met the criteria after a median of 3.4 months. Both models supported the same conclusion of a significant increase with time of volume loss in proximity to the ventricular interface.

Table 3

Parameter Estimates of the Nonlinear Mixed Effects Model Analyzing the Change in Thalamic Volume in MS Patients Compared to Children with Monophasic Demyelinating Disorders, as a Function of the Distance from the Ventricular Surface and Time Elapsed from Presentation

Parameter	Value	Standard Error	df	t	p
Intercept, I	1.67	1.43	980	1.16	0.25
Asymptote, a	2.60	0.80	980	3.26	0.0011
Slope, b	0.32	0.12	980	2.80	0.0052
Crossover point, c	8.18	0.57	980	14.31	<0.0001
Asymptote:timeC, a time	0.33	0.11	980	3.09	0.0021
Asymptote:timeC 2, a time2	0.01	0.05	980	0.15	0.88

Random effect: subject standard deviation = 4.88; residual standard deviation = 7.39; R 2 = 0.706. The significant interaction between time elapsed from clinical presentation and the asymptote (asymptote:timeC, a time) indicates a significant increase over time of the amount of volume change occurring in proximity to the cerebrospinal interface in MS compared to monophasic children.

MS = multiple sclerosis.

Table 4

Asymptote Estimate at Different Time Elapsed from Clinical Presentation

Years after Presentation	Asymptote	Standard Error	t	p	95% CI
1	3.01	0.800	3.77	0.0002	1.44–4.58
2	3.29	0.69	4.76	<0.0001	1.93–4.65
3	3.59	0.68	5.26	<0.0001	2.25–4.93
4	3.90	0.70	5.54	<0.0001	2.52–5.28

The asymptote estimate, capturing the difference between the ceiling and bottom values of the sigmoid model and therefore the amount of volume change occurring in proximity to the cerebrospinal fluid interface in multiple sclerosis compared to monophasic acquired demyelinating syndrome children, is already statistically significant at the earliest time point of 1 year from clinical onset. At the subsequent time points, it maintains its significance while increasing the magnitude of the observed effect.

CI = confidence interval.

To ensure that our analysis was not confounded by potential differences in the ability to accurately define the thalamic boundaries (ie, subtle differences between imaging characteristics of the CSF–thalamic interface versus the WM–thalamic interface), we repeated the above analyses after filling all nonthalamic voxels with the mean intensity of normal‐appearing WM. We observed distance‐related volume abnormalities in the MS and monoADS cohorts that were similar to those observed in the main analysis. The MS‐specific gradient at the CSF interface was again evident, even in the first year of clinical disease (p = 0.002).

Discussion

A surface‐in gradient of subpial cortical damage has been described in adult MS autopsy series and implicated as an important pathologic substrate of nonrelapsing progressive MS in late disease.5, 9, 38 A key question, however, is whether such a process also occurs subclinically, much earlier in the course of MS. The answer would have important implications for our fundamental concept of when treatments should be considered, as well as new treatment trials pursued, that target the nonrelapsing progressive MS biology.

The thalamus, which exhibits an early and disproportionate involvement in both adult onset and pediatric onset MS,17, 18, 28, 39 has 2 distinct interfaces: one with the WM and the other at the ventricular surface, where it may be exposed to a CSF‐mediated mechanism of injury. By using the resolution of voxel‐by‐voxel measurements of volume changes over time, and NLME modeling, we were able to compare the degree and trajectories of thalamic volume abnormalities in children with MS against the trajectories in controls, as a function of distance within the thalamus from the ventricular interface. Important features of our comparison were the inclusion of children with monophasic CNS demyelination, a population of patients that shares an inflammatory WM insult, but not the chronic biology of MS; and the inclusion of a healthy pediatric population, which enabled us to account for normal growth trajectories, as required for valid interpretation of volumetric data in the still‐developing pediatric brain.

Our study demonstrates for the first time the presence of a dual pattern of damage in the thalamus of children with MS. On one hand, we observed that atrophy of thalamic tissue adjacent to WM was a shared feature of both monophasic demyelination and MS.

This overlapping pattern of thalamic damage might reflect consequences of retrograde and/or anterograde degeneration of axons in WM pathways (such as spinothalamic and thalamocortical pathways) resulting from focal inflammatory lesions outside the thalamus, present in both children with MS and those with monoADS. An impact of such injury on normal thalamic growth trajectories of both children with MS and those with monophasic conditions is in keeping with our recent work showing that failure of expected maturation and progressive loss of WM integrity in the overall brain occurs not only in pediatric onset MS, but also in children with a single demyelinating attack involving the brain WM.23

Unique to the MS cohort in our study, we observed an additional pattern of thalamic atrophy, with greatest tissue losses in close proximity to the CSF interface and a decreasing magnitude of loss with distance from it. Although we focused our pediatric MS imaging study on the thalamus, we suspect that this graded injury pattern may manifest more generally in periventricular structures of MS patients. Consistent with this are recent imaging studies in adult onset MS reporting an abnormal and graded pattern of periventricular MTR abnormality,6, 7 with declining MTR levels observed in adjacent WM soon after the sentinel clinical event.8

In addition to our novel finding of a ventricular surface‐in gradient of thalamic injury that already exists at the time of initial clinical presentation of children with MS, our serial study further reveals ongoing volume losses in the same pattern at the ventricular–thalamic interface when these children were followed over time. Notably, a worsening gradient of thalamic abnormality was already discernable in the children with MS over the first year of follow‐up. Future work should focus on the relative rate of thalamic tissue loss (ventricular side vs WM side) through the generation of a dual model, able to capture the simultaneous changes associated with each of the two patterns of damage; on whether focal lesions within the WM and global or regional brain volume losses differentially correlate with each aspect of thalamic insult; and on whether the surface‐in trajectory of thalamic injury adjacent to the ventricular surface better correlates with measures of nonrelapsing progressive disease than does the trajectory of the thalamic injury on the WM side.

Our modeling (both parametric and semiparametric) of the trajectory of the pediatric MS‐specific surface‐in gradient of thalamic damage appears consistent with a diffusion‐mediated process emanating at the CSF interface.32 In such a process, the tissue in closest proximity to the CSF becomes saturated early and atrophies at a higher rate, and the effect exponentially decreases moving in depth, finally leveling off again where the concentration gradient no longer contributes to damage. The imaging evidence of such an injury gradient on the ventricular side of the thalamus of children with MS thus supports the theory put forth based on observations in adult MS that the subpial cortical damage may reflect the presence of one or more inflammatory, cytotoxic, and/or myelinotoxic soluble CSF factors, possibly released by immune cells within the CSF or situated within meningeal infiltrates.9, 40, 41, 42, 43 Several in vitro studies further support this hypothesis; cultured neurons exposed to CSF of MS patients (but not CSF from non‐MS controls) exhibited oxidative stress and axonal damage, with the implication of particular ceramides (C16:0 and C24:0) enriched in the CSF of the MS patients as possible mediators of this injury.13, 14 In other work, mixed CNS glial cells exposed to supernatants of B cells isolated from MS patients (but not to supernatants of B cells from HCs) exhibited increased death of both oligodendrocytes12 and neurons,16 and monoclonal recombinant antibodies derived from expanded B‐cell clones isolated from CSF of MS patients induced demyelination and astrocyte activation in spinal cord explants.44

Inherent to all studies of rare diseases, including pediatric onset MS, is the relatively small sample size. However, the magnitude of the CSF–thalamic effect we observed was sufficient to selectively implicate in the children with MS compared to the monoADS and HC groups. How imaging measures, including trajectories of thalamic injury, relate to clinical measures in such children will be the focus of future studies that will include larger numbers of participants appropriately powered for such assessments. In this regard, cognitive impairment is reported in approximately 30% of pediatric MS patients within the first 5 years of clinical disease onset,45 and we have previously shown a strong relationship between cognitive impairment and total thalamic volume.46 Of interest will be to correlate changes in cognitive performance with both the CSF‐adjacent and WM‐adjacent abnormalities of thalamic growth seen in pediatric MS. One might of course expect changes in thalamic volume (global, as well as CSF‐facing and WM‐facing) to account for only a portion of the variance in cognitive outcomes, and it will therefore be important to incorporate imaging measures of global and regional (supra‐ and infratentorial) volumes, cortical mantle thickness, and measures of WM integrity in such studies.

Although our study modeled thalamic gradients as a function of time since initial clinical presentation (and as we report in Table 1, the follow‐up time for the MS and monoADS groups was quite similar), one notes that there is no definitive way to align the exact dates of biological onset of those children with MS and the children with monoADS, as the biological onset of MS is now thought to precede (at least in the great majority of patients) their initial clinical presentation.17, 18, 47 A strength of the present study is that all patients were followed prospectively from initial clinical presentation with careful ascertainment of both clinical and MRI evidence of new disease activity, enabling reliable confirmation of when an initial clinical attack was the first clinical expression of MS. Although we cannot entirely exclude the risk of false‐negative diagnoses of monoADS (ie, individuals who may take a long time to have more disease activity to meet MS diagnostic criteria), this risk is likely to be minimal given the duration of follow‐up and the known natural history of pediatric MS, wherein the great majority of children with MS demonstrate clinical and/or MRI evidence of relapsing disease within 12 months of their first clinically evident attack.48, 49

Longitudinal measurements of volume change could be affected by several sources of bias,50 including biased registration, biased resampling space, and biased information transfer that relates to segmentation at multiple time points. Conventional registration is generally not symmetric and can introduce bias depending on which scan is chosen as the registration target, in the context of both linear51 and nonlinear52 registration. To minimize these possible biases in our study, we first created subject‐specific template spaces by linearly registering each scan to the MNI152 template for each participant; this avoids bias in the linear registration step and also provides an unbiased resampling space in which all further processing is performed. Then, each scan was registered to the baseline scan in the unbiased template space using a symmetric nonlinear registration algorithm (SyN; ANTs toolkit). It has been shown that use of symmetric registration, and the SyN algorithm in particular, limits biases introduced by conventional nonlinear registration, including in the multiscan serial case.52 Finally, we segmented the baseline scan, in the unbiased space, using an atlas‐based approach (ANIMAL; MINC toolkit). Our use of thalamic segmentation was only designed to provide a region of interest (not to make volume measurements), serving to limit possible biases related to segmentation at multiple time points.

A limitation of our study is that we do not assess the potential contribution of focal intrathalamic lesions, nor of potential differences between nuclei within the thalamus (whose different properties may make them more or less susceptible to injury), although a higher‐field magnetic resonance scanner would be required to more fully visualize such intrathalamic structures. The presence of nonuniform intensity within the thalamus allowed nonlinear registration to properly align local features, and thus provide local volume change measurements. Imperfections in the segmentation technique at the thalamic–CSF interface may have led to the inclusion of some nonthalamic (ie, CSF) voxels into the segmented thalamus. To minimize potential artifacts due to partial volume effects, we had excluded from the original analysis voxels in the thalamus that were within 2 mm of the CSF interface. One also notes that any erroneous inclusion of CSF voxels as thalamic voxels would have contributed to artifactual tissue “expansion” as the ventricular area is enlarging, and would therefore serve to cancel out, rather than enhance, the observation of preferential volume loss in proximity to the CSF boundary in patients with MS. Finally, our registration‐based volume change measurements use only a single segmentation per patient (designed to identify a region of interest) and hence are not sensitive to small errors in that segmentation, in contrast to segmentation‐based volume measurement techniques that measure volume directly on serial segmentations (and thus require highly accurate, consistent segmentation on each scan). The generation of more comprehensive modeling approaches, accounting simultaneously for the changes occurring at the two thalamic interfaces, would also allow assessment of the relative rates of change on each thalamic side in both the MS and monoADS cohorts simultaneously.

In summary, we demonstrated an evolving, MS‐specific, surface‐in pattern of thalamic injury that is present at the time of initial clinical presentation in children with MS. We postulate that this may reflect the very early presence of progressive (non‐relapsing) MS pathology in these children. Our study provides the rationale for considering that relapsing and non‐relapsing MS biological processes of injury overlap from very early in the disease course. Our modeling of the thalamic volume changes at the ventricular interface in children with MS at very early disease stages is compatible with the possibility that diffusible CSF factors mediate this injury, providing further impetus to study CNS‐compartmentalized processes and the precise cell and molecular events contributing to non‐relapsing MS damage. Additional work will be needed to determine whether the imaging and analytic approaches utilized here will prove useful as an outcome measure able to both discriminate and potentially monitor treatments that target progressive (non‐relapsing) MS mechanisms.

Author Contributions

G.F., R.A.B., and A.B.‐O. contributed to the conception and design of the study. G.F., R.A.B., R.M., B.A.‐B., E.A.Y., J.O., R.T.S., D.L.A., and D.L.C. contributed to the acquisition and analysis of data. G.F., R.A.B., A.B.‐O., B.B., R.A.M., and R.T.S. contributed to drafting the text and preparing the figures.

Potential Conflicts of Interest

Nothing to report.

Supplementary Table S1 Parameter estimates of the non‐linear mixed effects model analyzing the change in thalamic volume in MS patients compared to children with monophasic demyelinating disorders, as a function of the distance from the ventricular surface and time elapsed from MS diagnosis.

Click here for additional data file.