Multimodal magnetic resonance neuroimaging measures characteristic of early cART-treated pediatric HIV: A feature selection approach.

Children with perinatally acquired HIV (CPHIV) have poor cognitive outcomes despite early combination antiretroviral therapy (cART). While CPHIV-related brain alterations can be investigated separately using proton magnetic resonance spectroscopy (1 H-MRS), structural magnetic resonance imaging (sMRI), diffusion tensor imaging (DTI), and functional MRI (fMRI), a set of multimodal MRI measures characteristic of children on cART has not been previously identified. We used the embedded feature selection of a logistic elastic-net (EN) regularization to select neuroimaging measures that distinguish CPHIV from controls and measured their classification performance via the area under the receiver operating characteristic curve (AUC) using repeated cross validation. We also wished to establish whether combining MRI modalities improved the models. In single modality analysis, sMRI volumes performed best followed by DTI, whereas individual EN models on spectroscopic, gyrification, and cortical thickness measures showed no class discrimination capability. Adding DTI and 1 H-MRS in basal measures to sMRI volumes produced the highest classification performance validation accuracy = 85 % AUC = 0.80 . The best multimodal MRI set consisted of 22 DTI and sMRI volume features, which included reduced volumes of the bilateral globus pallidus and amygdala, as well as increased mean diffusivity (MD) and radial diffusivity (RD) in the right corticospinal tract in cART-treated CPHIV. Consistent with previous studies of CPHIV, select subcortical volumes obtained from sMRI provide reasonable discrimination between CPHIV and controls. This may give insight into neuroimaging measures that are relevant in understanding the effects of HIV on the brain, thereby providing a starting point for evaluating their link with cognitive performance in CPHIV.

INTRODUCTION

Without treatment, perinatal human immunodeficiency virus (PHIV) infection enters the central nervous system (CNS) soon after transmission. It is initially marked by elevated CSF white blood cell count and neopterin, followed by calcification of the basal ganglia, cortical and cerebral atrophy, and reactive gliosis in the maturing brain (George et al., 2009; Kieck & Andronikou, 2004; van Rie et al., 2007).

Combination antiretroviral therapy (cART) suppresses the elevated markers in each phase of infection and has improved clinical outcomes in children (Crowell et al., 2015; Cusini et al., 2013). However, the neurocognitive and behavioral profile of children on cART remains within the low average range (e.g., Garvie et al., 2014; Laughton et al., 2012; Nichols et al., 2016, 2015). A possible reason is differential crossing of the blood–brain barrier (BBB) by different antiretroviral (ARV) drugs (Caniglia et al., 2014; Letendre et al., 2008). Furthermore, the stage of brain development and progression of infection within which cART is initiated is crucial, as evidenced by the Children with HIV Early Antiretroviral (CHER) trial. This open label randomized controlled trial addressed when to start cART treatment in children with PHIV (CPHIV), and found that early time‐limited treatment (initiated before 12 weeks of age) reduced mortality, correlated with better neuropsychological test performance at 11 months, and provided better clinical outcomes than deferred continuous treatment (Cotton et al., 2013; Laughton et al., 2012; Violari et al., 2008).

A subset of children from the CHER trial were enrolled in a follow‐up longitudinal neuroimaging study with structural magnetic resonance imaging (sMRI), diffusion tensor imaging (DTI), and proton magnetic resonance spectroscopy (1H‐MRS) every 2 years from the age of 5 years. Compared with age and community‐matched uninfected children (controls), CPHIV at age 5 years showed various brain alterations, including larger subcortical gray matter structures and smaller corpora callosa (Randall et al., 2017), thicker cortex in bilateral frontal and left temporo‐insular regions, lower gyrification in left superior and bilateral medial orbitofrontal cortex (Nwosu et al., 2021), and regional fractional anisotropy (FA) decreases and mean diffusivity (MD) increases (Ackermann et al., 2016). CPHIV who started ART before age 12 weeks also had higher basal ganglia (bg) N‐acetyl‐aspartate (NAA) and total choline levels than controls (Mbugua et al., 2016). At age 7 years, among other findings, CPHIV had smaller volumes of gray and white matter (WM), bilateral hippocampus, putamen and right thalamus; lower gyrification in bilateral paracentral and right temporal regions, and thicker cortex in a small left inferior lateral occipital region (Nwosu et al., 2018); clusters with lower FA and higher MD in various WM tracts (Jankiewicz et al., 2017); but similar bg neurometabolite concentrations (Robertson et al., 2018) to controls.

These findings in the CHER children are consistent with other studies demonstrating HIV‐related irregularities in brain structure and WM integrity across childhood (Blokhuis et al., 2017, 2019; Hoare et al., 2019; Lewis‐de los Angeles et al., 2016, 2017; Wade et al., 2019; Yadav et al., 2017; Yu et al., 2019). Although a single MRI modality can identify pathological brain alterations, in isolation none can provide a comprehensive picture of the impact of HIV on brain structure and function. Many studies nowadays acquire data using multiple MRI modalities, but these data are rarely analyzed together to identify a signature of pediatric HIV. Due to the large number of measurements (features) generated in a typical neuroimaging study, it is a challenge to identify the most salient HIV‐related alterations.

Feature selection and supervised machine learning algorithms are increasingly applied to high‐dimensional neuroimaging datasets (e.g., Davatzikos, 2019; Janssen et al., 2018; Jollans et al., 2019; Woo et al., 2017; Wu et al., 2017) where a large number of neuroimaging features are combined and analyzed simultaneously. We aimed to build multimodal neuroimaging models, combining features from multiple MRI modalities, using elastic‐net (EN) regularization (Zou & Hastie, 2005) with repeated cross validation to identify MRI signatures of CPHIV on cART. A logistic EN regularization approach (Friedman et al., 2010; Zou & Hastie, 2005) uses penalized regression to perform variable selection, shrinking the weighting of unimportant variables to zero. Similar approaches were employed to classify dementia (Tohka et al., 2016) including Alzheimer's disease (Schouten et al., 2016) and Parkinson's disease (Bowman et al., 2016), as well as attention deficit hyperactivity disorder (Colby et al., 2012). However, no study to date has combined MRI modalities or used a multivariate pattern analysis approach to identify a set of neuroimaging features characteristic of CPHIV.

We hypothesized that a multivariate, multimodal approach combining neuroimaging data from sMRI, DTI and 1H‐MRS would discriminate better between 7‐year‐old controls and CPHIV on early cART than single modality analyses.

METHODS

Participants

Participants were 72 virally suppressed 7‐year‐old children from the CHER trial residing in Cape Town, South Africa (Cotton et al., 2013; Violari et al., 2008) who enrolled in a longitudinal follow‐on study and 55 age‐ and community‐matched uninfected control children. On the CHER trial, infants with PHIV who had CD4 percentage (CD4%) of at least 25% were randomized to one of three treatment arms. For one group, treatment was deferred until CD4% dropped below 25% in the first year or below 20% thereafter. The other two groups began cART immediately (before 12 weeks of age) with planned treatment interruption after, respectively, 40 and 96 weeks. Thereafter, cART was restarted if either CD4+ T‐cell declined to threshold levels or there was clinical evidence of severe HIV disease.

Due to evidence of accelerated disease progression following interruption, cART was not interrupted in a subset of the children who had initiated cART before 12 weeks as there was already evidence of disease progression. Fifteen infants with PHIV with CD4 < 25% were included in the early treatment arms. The cART regimen consisted of Zidovudine (ZDV) + Lamivudine (3TC) + Lopinavir‐ritonavir (LPV/r) (Cotton et al., 2013; Violari et al., 2008). Infants across all three groups showed first viral load suppression between 6 and 76 weeks of age. Sample characteristics of the children with PHIV are summarized in Table 1.

TABLE 1|

Sample characteristics of children with PHIV from the CHER trial.

CPHIV
Demographics
Observations (n)	72
Sex, Female, n (%)	36 (50%)
Age at scan (years)	7.20 (0.02)
Clinical measures at enrolment/baseline (6–8 weeks)
CD4 count (cells/mm3)	1800 (106.80)
CD4%	33 (1.21)
CD4/CD8	1.30 (0.08)
CD8 count (cells/mm3)	1682 (127.26)
CD8%	31 (1.24)
Viral load (copies/ml), n (%)
High (>750 000)	43 (58%)
Low (400–750 000)	29 (40%)
Suppressed (<400)	0
Clinical measures at scan
CD4 count (cells/mm3)	11652 (55.16)
CD4%	37 (0.73)
Viral load (copies/ml), n (%)
High (>750 000)	0
Low (400—750 000)	5 (7%)
Suppressed (<400)	67 (93%)
Treatment‐related measures
cART initiation before 12 weeks, n (%)	54 (75%)
Children with cART interruption, n (%)	40 (55.56%)
Age at cART interruption (weeks) a	70.40 (4.37)
Duration of cART interruption (weeks) b	61.44 (13.43)

Note: Values presented are “mean (SE)”, unless otherwise stated.

Based only on the children in whom treatment was interrupted; median 51.27 weeks, interquartile range 55.43 weeks.

Based only on the children in whom treatment was interrupted; median 32.71 weeks, interquartile range 33.07 weeks.

Since the focus of the current study is to identify neuroimaging features that distinguish CPHIV (irrespective of treatment arm) from children without HIV, we do not present a comparison of clinical measures for the three CHER treatment arms. Moreover, clinical measures of the three treatment groups have been presented in previous papers (e.g., Laughton et al., 2018; van Wyhe et al., 2021). Briefly, among the CPHIV in the current sample, enrolment (age 6–12 weeks) CD4, CD8 and viral load were similar across the three treatment groups (all p's > 0.06; pairwise Wilcoxon rank sum tests), as were CD4 and CD8 around age 7 years (p's > 0.10). In the early treatment arm with planned interruption at 96 weeks, five children had detectable viral loads around age 7 years, while all children in the other two arms were virally suppressed (<400 RNA copies/ml).

Image acquisition

Children were scanned around age 7 years on a 3 Tesla Allegra MRI scanner (Siemens Erlangen, Germany) at the Cape Universities Brain Imaging Centre (CUBIC) using a single channel head coil. Neuroimaging was performed in a single session lasting under 60 min without sedation in accordance with protocols approved by the Human Research Ethics committees of the Universities of Cape Town and Stellenbosch. Parents or guardians provided written informed consent and oral assent was provided by the children, who were first familiarized with the procedure in a mock scanner.

Structural and diffusion tensor imaging were performed as well as single‐voxel 1H‐MRS. The high‐resolution structural ‐weighted acquisition used a 3‐dimensional echo planar imaging (EPI) navigated multi‐echo magnetization prepared rapid gradient echo (MEMPRAGE) sequence (van der Kouwe et al., 2008) with field of view (FOV) , voxel size , flip angle , echo times (TEs) , repetition time (TR) , inversion time (TI) , bandwidth , matrix size 224 × 168, and 144 sagittal slices. Two diffusion‐weighted (DW) datasets were acquired with opposite phase encodings (anterior–posterior and vice versa) for EPI distortion correction during processing using a volumetric navigated (Alhamud et al., 2012) twice refocused spin‐echo (TRSE) sequence: , , 72 slices, matrix size 112 × 112, TR = 10,000 ms, , , , 30 diffusion directions. 1H‐MRS were obtained from 3 voxels of interest in the basal ganglia (bg), midfrontal gray matter (mfgm) and peritrigonal white matter (pwm) using an EPI navigated point resolved spectroscopy (PRESS) sequence (Hess et al., 2011): and chemical shift selective (CHESS) (Haase et al., 1985) water suppression; an acquisition without water suppression was performed in each voxel for water referencing.

Image pre‐processing and feature creation

For the structural T1‐weighted dataset, FreeSurfer version 6.0 (Dale et al., 1999; Fischl et al., 2002) was used for automated cortical reconstruction and volumetric segmentation. Skull‐strip corrections (removal of extra‐cerebral voxels) were performed if required as described by Nwosu et al. (2018). Thirty eight regional volumes from FreeSurfer's automatic segmentation were included as candidate features in our sMRI volumes feature set (Figure S1), and cortical thickness and local gyrification index (LGI) for each of the 34 regions of the Desikan–Killiany (DK) atlas (Desikan et al., 2006) in each cerebral hemisphere (LH & RH) for the sMRI cortical thickness and LGI feature sets, respectively (Figure S2). LGI's were included as a candidate feature set as LGI alterations have previously been reported in CPHIV (Lewis‐de los Angeles et al., 2016; Hoare et al., 2019; Nwosu et al., 2018).

For the DW dataset, field distortion, motion correction, and a voxel‐wise calculation of the diffusion tensor were done in Tortoise v.2.5.2. Voxel‐wise diffusivity measures were calculated in AFNI (Cox, 1996) as described in Jankiewicz et al. (2017). Average values of radial diffusivity (RD), axial diffusivity (AD), mean diffusivity (MD) and fractional anisotropy (FA) were calculated for each of the 20 white matter tracts in the Johns Hopkins University (JHU) atlas (Mori et al., 2005; Figure S3).

For spectroscopy, processing was performed as described by Robertson et al. (2018). Concentrations of 14 metabolites and their ratios to creatine including N‐acetylaspartate (NAA), creatine (Cr), choline (Cho), myo‐inositol (MI), lipids and macromolecules (MM) (Table S1) were obtained for each of the three SVS voxels using LCModel version 6.1 (Provencher, 2001).

Quality‐control measures were applied to each of the seven MRI‐derived feature sets: sMRI volumes (1st set), sMRI cortical thickness (2nd set), sMRI gyrification (3rd set), DTI (4th set), 1H‐MRS bg (5th set), 1H‐MRS mfgm (6th set), and 1H‐MRS pwm (7th set). Extreme outliers on measures from the T1‐weighted structural data and DW images with large motion artefacts and dropout slices were eliminated from analysis. Datasets with failed alignment to the JHU and DKT templates in standard space were excluded. Similarly, poor quality spectra were eliminated based on low SNR (SNR <6) or broad linewidth (FWHM>0.07 ppm) as calculated in LCModel version 6.1. Because some subjects did not complete the full scanning session and in some subjects not all modalities met quality control criteria, all seven feature sets were not present for all subjects. Subjects were therefore excluded from analyses for which they did not have the required feature sets.

Statistical analyses

We chose to use the embedded feature selection of EN logistic regression to identify a collection of neuroimaging features characteristic of pediatric HIV infection. We did this initially within each neuroimaging feature set separately (namely sMRI volumes, sMRI cortical thickness, sMRI gyrification, DTI, 1H‐MRS bg, 1H‐MRS mfgm, and 1H‐MRS pwm) before combining features from different sets. EN models were implemented in the R programming language version 3.5.3 (R Core Team, 2018).

Confounders that may affect imaging measures may be included as predictors in a predictive model, that is, estimated jointly with neuroimaging data (Linn et al., 2016) or the effects of possible confounders may be regressed out of the neuroimaging data before the model is constructed (e.g., Adeli et al., 2018). Choosing the former approach, we included sex and age at scan as confounders for all EN models with no shrinkage penalty, and total intracranial volume (TIV) as an added confounder for models that included sMRI volumes.

Classification performance

To estimate how accurately each classification model performs on an independent data set, we implemented k‐fold cross‐validation (CV) (Kohavi, 1995). We chose to use k = 10 folds as is commonly used in neuroimaging studies which have relatively small n since it avoids overfitting and selection bias (Colby et al., 2012; Schouten et al., 2016; Tohka et al., 2016).

Initially, for sMRI volumes, sMRI cortical thickness, sMRI gyrification, DTI, 1H‐MRS bg, 1H‐MRS mfgm, and 1H‐MRS pwm separately, the classification performance was assessed by computing accuracy, specificity, and sensitivity at the optimum cut‐off point, as well as the mean area under the receiver operator characteristic curve (AUC) across 10 CV folds (Figure 1) and 100 iterations (Krstajic et al., 2014; Schouten et al., 2016). We identified the best performing single modality model as the one with highest AUC.

FIGURE 1

Illustration of classification and feature selection using repeated 10‐fold cross validation employed in this study.

Feature identification

To identify features that typify cART‐treated pediatric HIV infection and to assess the stability of the selected features, the EN model that gave the highest AUC was refitted on the whole sample (i.e., all subjects) using optimized parameters obtained during the 100 training iterations. The relative importance of a selected feature was evaluated in two ways: the frequency (f) of selection with a non‐zero weighting () across the 100 repeats, and the average value of its weighting () in the final model over all iterations. Features were determined to be relevant if their frequency of occurrence was >75% (selected at least 75 out of 100 retraining iterations). This was an empirically determined threshold that produced the most overlap in features between the single modality and multimodal models. The same features were selected when the models were trained once on the entire data set.

Construction of multimodal models

To create multimodal models, the entire set of features from the best performing single modality model (highest AUC) was concatenated pairwise with each of the remaining six feature sets and more EN‐penalized regressions were run. To allow unbiased selection of multimodal features, at each stage we included all the features in each set. We then computed the classification performance and the combination resulting in the largest AUC was concatenated in turn with each of the remaining feature sets. This process was repeated until there was no further increase in AUC from concatenating the features of the best performing combination with another feature set.

To compare the best performing multimodal and single modality models we used the Akaike information criterion (AIC) (Akaike, 1973; Anderson & Burnham, 2002). This allowed us to evaluate which model provided the best trade‐off between simplicity (smaller number of parameters) and goodness of fit (minimizing information loss). We then chose the multimodal combination with minimum AIC. For large datasets, minimizing AIC is equivalent to minimizing the CV error (Peng et al., 2013; Stone, 1977).

We then identified multimodal features for the combination resulting in the minimum AIC using the same procedure described for the individual modality models in Section 2.4.2.

RESULTS

Number of observations used for analyses

Sample characteristics of participants whose data survived quality control for each modality are listed in Table 2. Structural imaging feature sets had the largest number of observations (n = 125), and spectroscopy in the basal ganglia had the fewest (n = 79). The two CPHIV excluded from the sMRI feature sets only had usable DTI measures. The number of observations for multimodal feature sets was necessarily smaller (lower number of subjects n) than the largest individual feature sets from structural imaging.

TABLE 2

Sample characteristics of participants (n) in each modality (and combination of modalities) after quality control.

Feature set		CPHIV	Controls
sMRI measures (n = 125) a  Volumes Gyrification Thickness	Observations	70	55
	Females, n (%)	35 (50%)	24 (44%)
	Age at scan (se)	7.20 (0.02)	7.24 (0.02)
1H‐MRS mfgm (n = 105) b	Observations	60	45
	Females, n (%)	31 (52%)	19 (42%)
	Age at scan (se)	7.19 (0.01)	7.24 (0.02)
DTI (n = 104) c	Observations	59	45
	Females, n (%)	31 (53%)	21 (47%)
	Mean age (se)	7.22 (0.02)	7.24 (0.02)
1H‐MRS pwm (n = 99) d	Observations	58	41
	Females, n (%)	30 (52%)	18 (44%)
	Age at scan (se)	7.19 (0.01)	7.24 (0.02)
1H‐MRS bg (n = 79) e	Observations	47	32
	Females, n (%)	23 (49%)	15 (47%)
	Age at scan (se)	7.23 (0.01)	7.19 (0.02)
sMRI volumes + DTI (n = 102) f	Observations	57	45
	Females, n (%)	30 (53%)	21 (47%)
	Age at scan (se)	7.21 (0.02)	7.24 (0.02)
sMRI volumes + DTI + 1H‐MRS bg	Observations	40	29
	Females, n (%)	24 (60%)	16 (55%)
(n = 69) g	Age at scan (se)	7.21 (0.02)	7.24 (0.03)

Twenty‐three subjects had no DTI measures, 46 had no 1H‐MRS bg, 21 had no 1H‐MRS mfgm, and 27 had no 1H‐MRS pwm.

Fourteen subjects had no DTI measures, one did not have structural measures, 31 had no 1H‐MRS bg, and 14 had no 1H‐MRS pwm.

Thirty‐five subjects had no 1H‐MRS bg, 13 had no 1H‐MRS mfgm, 17 had no 1H‐MRS pwm, and two had no structural measures.

One had no structural imaging, 12 had no DTI, 29 had no 1H‐MRS bg, and eight had no 1H‐MRS mfgm.

All had structural measures, 10 had no DTI, five had 1H‐MRS mfgm, and nine had no 1H‐MRS pwm.

All had structural and DTI, 33 had no 1H‐MRS bg, 1H‐MRS bg, 12 had no 1H‐MRS mfgm, and 16 had no 1H‐MRS pwm.

All had structural measures, DTI, and 1H‐MRS bg, five had no 1H‐MRS mfgm, and six had no 1H‐MRS pwm.

Classification performance for individual feature sets

Building EN models using data from individual modality feature sets, we found that structural volumes and DTI metrics produced better discrimination between CPHIV and controls than other measures considered (summarized in Table 3).

TABLE 3

Classification performance evaluation metrics of the seven MRI derived feature sets.

Modality	n	AUC	Sensitivity (%)	Specificity (%)	Accuracy (%)
sMRI volumes	125	0.71	81	69	74
sMRI thickness	125	0.55	80	60	71
sMRI gyrification	125	0.54	75	65	65
1H‐MRS mfgm	105	0.58	77	60	61
DTI all regions	104	0.62	75	69	72
1H‐MRS pwm	99	0.49	72	32	59
1H‐MRS bg	79	0.58	64	69	61

Note: We present the AUC (degree of separability) along with sensitivity, specificity, and accuracy for each individual classifier. The number of observations (n) after quality control is also given.

Selected features for the best performing single modality

We show the 12 relevant features from the sMRI volume model in Figure 2 and the ranking of these features in Table 4. The sign of the coefficients indicates smaller volumes of the relevant structures except for CSF in children living with HIV than in controls. Figure 3 shows features selected in the DTI model and their feature ranking in Table S2.

FIGURE 2

Anatomical locations of the relevant features from the sMRI volume feature set. The sMRI volumes were from an automatic segmentation with FreeSurfer (Fischl et al., 2002). CWM, cerebellar white matter; CC, corpus callosum; CSF, cerebrospinal fluid).

TABLE 4|

Relevant features of the sMRI volumes ranked by the absolute value of their average weighting across 10 folds.

CV penalized EN regularization model
Relevant features	Frequency	β¯ × 10−2 a
Right vessel	100	−2.37
Optic Chiasm	98	−0.64
Left Amygdala	99	−0.25
Left Pallidum	100	−0.20
Right Pallidum	100	−0.18
Right choroid plexus	95	−0.15
Right amygdala	93	−0.08
Posterior corpus callosum (CC)	78	−0.07
Cerebrospinal fluid (CSF)	93	0.05
Right ventral diencephalon (DC)	84	−0.05
Left ventral DC	77	−0.02
Right cerebellum white matter (CWM)	90	−0.01

Mean weighting of the feature in the penalized EN regularization model. A negative mean weighting indicates a smaller measure in CPHIV than in controls, and adjusted for age at scan, sex, and TIV.

FIGURE 3

Selected relevant DTI tracts. The regions of interest for the DTI feature set were from the John Hopkins University (JHU) atlas (Mori et al., 2005).

Classification performance for multimodal feature sets

The stepwise concatenation procedure that began with sMRI volumes (AUC = 0.71) gave a multimodal feature set consisting of sMRI volumes + DTI + 1H‐MRS bg with AUC = 0.80 and 85% accuracy (Table 5). Table S3 provides accuracy, sensitivity, and specificity for the best combination in each step. No improvement of AUC was seen with further addition of feature sets.

TABLE 5|

Multimodal classification performance in steps of improving AUC (in bold).

Concatenation steps	sMRI volumes	DTI all regions	1H‐MRS bg	1H‐MRS mfgm	sMRI thickness	sMRI LGI	1H‐MRS pwm
1. —	0.71	0.62	0.58	0.58	0.55	0.54	0.49
2. sMRI volumes (n = 125)	—	0.78	0.75	0.68	0.61	0.60	0.62
3. sMRI volumes + DTI (n = 102)	—	—	0.80	0.78	0.72	0.65	0.76
4. sMRI volumes + DTI + 1H‐MRS bg (n = 69)	—	—	—	0.77	0.73	0.64	0.76

Note: The first column gives the feature set combination with which the feature set in each column is concatenated. Only subjects with data for all measurements were included. Since the number of subjects n was different for each multimodal feature set, each concatenation step was performed on a smaller sample than the previous one.

Goodness of fit comparison for single vs multimodal models

We compared our models in terms of predictive ability (CV) and fit to the data (AIC) and found that although the model with the lowest CV error was sMRI volumes + DTI + 1H‐MRS bg () in Table 5, the simpler multimodal model comprising only the sMRI volumes and DTI feature sets had a lower AIC (; Table 6)—making it the simplest model to minimize information loss.

TABLE 6

Comparison of non‐nested multimodal classifiers with AIC.

Candidate models	k a	AIC_c	Relative likelihood e^AIC_cmin−AIC2
sMRI volumes	16	99.60	0.013
DTI	16	102.87	0.002
1H‐MRS bg	29	163.12	2.056 × 10−16
sMRI volumes + DTI	26	90.88	1.000
sMRI volumes + DTI + 1H‐MRS bg	28	102.62	0.002

Note: Second order AIC (AICc) because ratio of sample size to number of estimated parameters was <40.

Number of estimated parameters in the model = relevant features + confounders + deviance (measure of goodness of fit).

To confirm the sensitivity of the model to distinguish virally suppressed CPHIV from controls, we repeated the analyses excluding the five children with detectable viral loads. Similar classification metrics (sensitivity, specificity, accuracy, and AUC) were obtained. Further, repeating the analyses of the three best performing feature sets (sMRI volumes, DTI, and sMRI volumes + DTI) using only the subjects common to all feature sets (n = 102) yielded similar classification performance measures (Table S4), confirming the stability of the models.

Finally, we examined for the three best performing feature sets whether nadir CD4 and nadir CD4%, which reflect historical processes and are known to be strong indicators of brain integrity, and differences in treatment interruption within our sample, biased model performance. We found no differences in nadir CD4, nadir CD4% (Wilcoxon Rank Sum tests, all p's > .22), nor the percentage of children with treatment interruption (Chi‐square test p's > .53), between correctly and incorrectly classified CPHIV.

Features selected in the multimodal feature set

We found 22 stable features from the multimodal combination of volumes and DTI (ranked in Table S5). Although most selected predictors in the multimodal combination were the combined predictors selected in the two separate single modality EN models, seven less relevant features, namely right choroid plexus, posterior corpus callosum (CC), right ventral diencephalon (DC), MD, RD and AD in left corticospinal tract, and RD in left cingulum, were replaced by four features in the multimodal combination that were not selected in the single modality analyses: volumes of brainstem and right lateral ventricle, RD in right cingulum, and AD in the temporal part of the right superior longitudinal fasciculus (illustrated in Figure 4).

FIGURE 4|

Four extra features that appear in the multimodal analysis and were not part of the single modality analyses of sMRI volumes and DTI. These are in addition to all except three (right choroid plexus and ventral diencephalon, and posterior corpus callosum) of the sMRI features in Figure 2 and the DTI tracts, except left corticospinal tract, in Figure 3. The DTI regions of interest are from the John Hopkins University (JHU) atlas (Mori et al., 2005) and the sMRI volumes from FreeSurfer (Fischl et al., 2002).

DISCUSSION

Feature selection and classification performance

In line with our hypothesis, multimodal neuroimaging improves classification performance compared with single modalities. We identified a neuroimaging signature of cART‐treated pediatric HIV comprising 22 features from a multimodal combination of regional volumes and DTI measures in several tracts—these include smaller bilateral globus pallidus and amygdala, as well as increased radial and mean diffusivities in the uncinate fasciculus and right corticospinal tract. We infer that there are more consistent effects of cART‐treated HIV on subcortical structures and WM microstructure at age 7 years for these children than on cortical thickness, gyrification, and neurometabolite concentrations in the basal ganglia, midfrontal gray matter, and peritrigonal white matter.

Combining sMRI and DTI previously improved classification performance in other pathologies like Alzheimer's disease (Schouten et al., 2016), schizophrenia (Sui et al., 2015), depression (Kambeitz et al., 2017), and Parkinson's disease (Bowman et al., 2016). While we and these authors relied primarily on AUC as an overall measure of model performance, other studies in dementia (Tohka et al., 2016), mild cognitive impairment (Ritter et al., 2015), attention deficit hyperactivity disorder (Colby et al., 2012), and schizophrenia (Cabral et al., 2016) reported increases of 2–7% in classification accuracy with a multimodal approach. Similarly, our initial classification accuracies of 74%, 72%, and 61% for volumes, diffusion measures, and basal ganglia metabolites, respectively, improved to 85% for the multimodal combination. Previous studies could discriminate adults living with HIV from controls with accuracies >80% using classifiers incorporating embedded feature selection, namely LASSO on DTI measures (Tang et al., 2017), and chained regularization on sMRI (Adeli et al., 2018). The combination of multimodal MRI measures with clinical features has also been found to predict neurocognitive impairment in adults with HIV better than clinical features or MRI measures alone, suggesting that these techniques have the potential to be clinically useful (Xu et al., 2021).

This is the first study to combine MRI modalities and use a multivariate pattern recognition approach to identify features characteristic of pediatric HIV infection. This may give insight into neuroimaging measures that are relevant in understanding the effects of HIV on the brain, thereby providing a starting point for evaluating the link with cognitive performance in CPHIV. Although most of the predictors of cART‐treated pediatric HIV in the multimodal combination are simply the combined relevant features selected in each single modality when analyzed separately, the multimodal EN model does not retain redundant features (Chandrashekar and Sahin, 2014; Khaire and Dhanalksmi, 2019; Mwangi et al., 2014). We found the volumes of the posterior corpus callosum, right choroid plexus and right ventral diencephalon, as well as MD, RD and AD in left corticospinal tract, and RD in right cingulum to be relevant in the single modality analyses, but redundant in the multimodal analysis. Similarly, brainstem and right lateral ventricle volumes, RD in right cingulum, and AD in the temporal part of the right superior longitudinal fasciculus were relevant only in the multimodal analysis. This implies that when sMRI volumes and DTI are combined, the latter four features have better, or at least equivalent, explanatory power compared with the seven dropped single modality features, avoiding model overfitting (cross validation), and minimizing information loss.

We used both cross‐validation and the second order AIC (Anderson & Burnham, 2002) to evaluate our models in terms of predictive ability and fit to the data. Although these are equivalent in certain cases (Peng et al., 2013; Stone, 1977), in our model parameters were optimized to minimize cross‐validation error rather than AIC. Interestingly, the complex model (sMRI volumes + DTI + 1H‐MRS bg) had the lowest CV error while the simpler model (sMRI+ DTI) provided the best compromise between simplicity and minimizing information loss, that is, minimum AIC. Since our data set was relatively small for the complex model (n = 69), our training and tests were also necessarily small, which may reduce reliability of cross‐validation (Varoquaux, 2018). It is worth testing whether the classification performance from this study would be similar if the models were used to predict HIV status in children from other cohorts. If our data are similar to other CPHIV cohorts, we would expect our results to remain valid using the model that minimizes CV error. However, if the distribution of our data is different to the actual data distribution of other cohorts, the models with fewer parameters (and lower AIC) might perform better.

Pediatric neuroimaging research is frequently limited by number of observations (e.g., Benki‐Nugent and Boivin et al., 2019; George et al., 2015; van Arnhem et al., 2013; Wade et al., 2019), since relatively small cohorts of subjects are scanned. Our study sample size of 127 was further reduced after quality control criteria had been applied to each modality. However, larger sample size does not directly translate to higher classification performance, since models using DTI features (n = 104) showed better performance than those using cortical thickness and gyrification (largest n = 125); furthermore, multimodal classifiers (n = 102 and n = 69) have better performance that those from any single modality. Therefore, we infer that differences in sample size between models did not play a substantial role in determining classification performance. This was confirmed by showing that the classification performance of the single modality sMRI volumes and DTI analyses remained essentially unchanged when repeated with only the subjects (n = 102) common to the sMRI volumes + DTI model (Table S4).

In addition, as there were more CPHIV than controls, prediction models are biased toward higher sensitivity than specificity. Cross validation and random shuffling of subjects during repeat cycles were performed to minimize bias, however higher sensitivity than specificity suggests that these models still tended to identify CPHIV better than controls. Ninety‐three percent of the CPHIV had suppressed plasma HIV RNA; exclusion of the five children with detectable viral loads did not lower classification performance.

There is still debate about how to incorporate confounders into predictive models, depending on whether the goal of the model is primarily biomarker development, disease diagnosis or to understand the mechanisms and neuroimaging patterns that characterize the disease (Linn et al., 2016). Applying feature selection to all variables including possible confounders such as age, sex and total intracranial volume with a shrinkage penalty showed that none of these features predicted HIV status. We conclude that these are not confounders in this sample, but since these variables may be confounders in the larger CPHIV population, they were included in the model with no shrinkage penalty. Therefore, other coefficients in the models are adjusted for these variables. Furthermore, as we found no differences in nadir CD4 count, nadir CD4%, nor the percentage of children with treatment interruption, between correctly and incorrectly classified CPHIV, we conclude that these historical factors were not driving the findings.

Including non‐imaging confounders—in this case sex and age, is possibly sub‐optimal for small sample sizes in which there may be a bias in the training sample—a relationship between a confounder and the outcome variable that is not representative of the wider population (Linn et al., 2016). In our overall data set, there is no difference in age distribution between CPHIV and control children. However, there are on average 7% more females in the CPHIV groups (Table 2). If non‐imaging variables are associated with the imaging measures in a disease‐dependent way, estimating disease patterns jointly with non‐imaging variables may reveal the relative importance of neuroimaging versus demographic and clinical features in the classification (Linn et al., 2016). However, this can be a problem when the distribution of confounders in the training sample does not match the distribution in the population of interest (Linn et al., 2016; Rao et al., 2017).

Implications of selected features in cART‐treated pediatric HIV

Previous studies in this cohort have also reported differences in subcortical volumes in cART‐treated pediatric HIV (Nwosu et al., 2018; Randall et al., 2017). However, with training and cross‐validated feature selection, a feature that alone is not highly predictive of HIV can be discovered when combined with another feature or set of features. This may explain measures in our EN regularization models that were not altered in previous investigations in the same cohort, for example right lateral ventricle volume. Furthermore, the EN analysis allows us to consider many candidate features previously unexamined because of the need to limit multiple comparisons. For example, we found reduced volume of bilateral amygdala and expansion of CSF including third and fourth ventricles, structures that were not included in previous analyses in this cohort (Nwosu et al., 2018; Randall et al., 2017).

Although Nwosu et al. (2018) found reduced putamen volumes in the CHER cohort at age 7 years, the putamen is not a relevant feature in our single modality or multimodal EN models in the same population. This suggests that smaller putamen may not be a robust finding across the CPHIV population and may reflect a change that is specific to a CPHIV subgroup, for example, children who started cART after 12 weeks. In fact, Blokhuis et al. (2017) reported trend‐level putamen increases in CPHIV children, and Randall et al. (2017) similarly reported larger putamen volumes at age 5 in the CHER cohort that were associated with delayed cART initiation.

Our results suggest that HIV may affect globus pallidus more than other components of the basal ganglia. In a cohort of Thai CPHIV, there was shape deformation of the right pallidal surface and baseline CD4 count was negatively associated with left globus pallidus volume (Wade et al., 2019). Shape deformation has also been negatively associated with peak viral load in adolescents (Lewis‐de los Angeles et al., 2016). Although our models suggest reduced volume of the left and right globus pallidus in CPHIV at this age, larger globus pallidus volumes are reported for CPHIV aged 5 in the CHER cohort (Randall et al., 2017) and CPHIV <12 years, but not in older children (Paul et al., 2018). The effect of HIV on the globus pallidus may therefore vary with factors such as age and past HIV disease severity.

Ventricular enlargement is among the more common neuroimaging findings in children with HIV on cART (van Arnhem et al., 2013). Although total CSF has not previously been examined in this cohort, our finding of increased CSF volume (collective expansion of the ventricles, that is, lateral ventricles + third ventricle + fourth ventricle) in CPHIV is therefore not unexpected. Surprisingly, despite increased total CSF, our EN models suggest reduced volume of the right lateral ventricle in HIV infection, despite previous analysis of this data showing no difference in volumes of the lateral ventricles in CPHIV (Nwosu et al., 2018). Previous studies have shown that children living with HIV have lateral ventricles that are asymmetrical (Bruck et al., 2001; Thompson and Jahanshad, 2015) and dilated, with a reduction in total brain volume (Mariam & Assefa, 2012).

Previous studies in cART‐treated children found increased MD in superior longitudinal fasciculus at ages 8–12 years (Hoare et al., 2012) and 13 to 17 years (Li et al., 2015), while Uban et al. (2015) found increased MD in the inferior longitudinal fasciculus between 16 and 20 years of age. Others described decreased AD in left superior longitudinal fasciculus in cART‐treated children aged 6–15 years, associated with poorer fronto‐striatal cognition (Hoare et al., 2015). We find only increased AD in the right temporal part of the superior longitudinal fasciculus. MD and RD differences along the cingulum have been reported in adults living with HIV (Zhu et al., 2013), however we were unable to find studies that report on the cingulum in cART‐treated children.

We noted increased RD and MD in the right corticospinal tract, suggestive of reduced myelination (Lebel & Beaulieu, 2011) and less organized structure (Feldman et al., 2010), respectively. At age 5, Ackerman et al. (2016) found two clusters in right corticospinal tract with lower FA predominantly attributable to increased RD. Increased MD in the same tract was due to both increases in RD and AD. In a follow‐up study at age 7 years, Jankiewicz et al. (2017) found no continued WM damage in this tract, possibly attributable to cART. Similarly, we found increased RD and MD in bilateral uncinate fasciculus, and increased AD in right uncinate fasciculus, as relevant features of cART‐treated HIV, although no voxel‐wise group differences were found in any clusters that overlap these tracts at age 7 years (Jankiewicz et al., 2017).

The inconsistency between these results and ours may be explained by noting that in this study we employed a region‐of‐interest (ROI) method instead of voxel‐based group analyses as employed by Jankiewicz et al. (2017) and Ackerman et al. (Ackermann et al., 2016). Where regions affected by HIV are substantially smaller than the atlas regions, differences found in voxel‐wise analysis might not be detectable when averaging over an ROI. Conversely, averaging across voxels in a noisy ROI could also provide more power to detect differences not apparent in a voxel‐wise comparison. This may also explain why, although Nwosu et al. (2018) and others found reduced cortical gyrification in HIV (Lewis‐de los Angeles et al., 2016; Hoare et al., 2019), gyrification was not an important feature in our models.

Finally, spectroscopy measures in the basal ganglia, midfrontal gray matter, and peritrigonal white matter are poor predictors of HIV at age 7 years according to our individual modality EN analyses. Correspondingly, in the same study population, Robertson et al. (2018) using linear regression models found no difference in metabolite concentrations in the basal ganglia at age 7. In our study, data were analyzed cross‐sectionally at age 7, which prevented us from exploring longitudinal changes due to cART‐treated HIV. These children were followed and scanned at ages 5, 7, 9, and 11 years, and are thus well characterized. We recommend a future longitudinal multimodal neuroimaging study. Lifelong exposure to HIV and cART in this population warrants a longitudinal assessment of these signatures, their changes over time, relationship to neuroinflammation, neurocognitive impairment, and vascular dysfunction and utility in predicting neurocognitive outcomes.

Other relevant predictors of cART‐treated HIV

In addition to the neuroimaging features described above, some rarely examined features that are specific to FreeSurfer segmentation were found to be relevant predictors in our models. These include right vessel of the basal ganglia, left ventral diencephalon, optic chiasm, and the brainstem.

The left and right vessel are small labels in FreeSurfer for vessels entering the basal ganglia that would otherwise be labelled as lesions (Filipek et al., 1994; Fischl et al., 2002). Reduced volume of the right vessel in CPHIV may reflect altered blood flow to the basal ganglia, which has previously been detected in cART‐treated perinatally HIV‐infected children with mild psychomotor symptoms (Blokhuis et al., 2017).

The FreeSurfer ventral diencephalon label contains the hypothalamus, its main component, along with basal forebrain, sub‐lenticular extended amygdala, and a large portion of the ventral tegmentum (Makris et al., 2008; Makris et al., 2013). As volume of the hypothalamus is difficult to measure in vivo (Makris et al., 2013), other investigators have used this FreeSurfer‐defined ROI as the best representation of the neuroendocrine or homeostatic system in pathologies such as pediatric neurotrauma (Bigler et al., 2019), major depression disorder (Kim et al., 2019), and first‐episode schizophrenia (Emsley et al., 2015). Endocrine dysfunction is extensively reported in children living with HIV (Chantry et al., 2007; Loomba‐Albrecht et al., 2014; Pandey et al., 2018; Rondanelli et al., 2002), however, no previous neuroimaging studies measured volume of hypothalamus or other parts of the ventral diencephalon in people living with HIV.

The optic chiasm is a small visual pathway—part of the brain where optic nerves partially cross. It is located at the bottom of the brain immediately inferior to the hypothalamus, and is essential for binocular vision and depth perception (Herrera et al., 2019). Although in the pre‐cART era between 40 and 90% of adults living with HIV had ocular involvement (Govender et al., 2011), the prevalence of ocular manifestations of HIV infection in children is very low (Dennehy et al., 1989; Padhani et al., 2000) suggesting that this finding is unlikely to have any functional consequence as is the case in other studies on optic chiasm abnormalities. Brainstem morphometry is not often examined. The brainstem includes the sensory pathway for the auditory nerve. Many studies have reported hearing impairment in children living with HIV on cART (Chao et al., 2012; Hrapcak et al., 2016; Torre et al., 2012; Torre III, 2016), which are commonly associated with abnormalities in auditory brainstem responses. The first relay of the primary auditory pathway occurs in the cochlear nuclei in the brainstem, where the basic decoding of duration, intensity, and frequency occurs.

Some of these selected features may not on their own be highly predictive of HIV infection, but their specificity to HIV may be revealed in combination with another feature or set of features. Discovery of common patterns of features characteristic of HIV has the potential to reveal mechanisms underlying the effects of HIV on the developing brain. Although it is unclear how the selected features in this study are related to each other, some seem to have anatomical links. For example, the uncinate fasciculus connects parts of the limbic system, which includes amygdala and diencephalon, with the orbitofrontal cortex. The cingulum also allows components of the limbic system to communicate, and the large fiber bundles of the corticospinal tract connect the cortex with spinal cord via the brainstem. The selection of these features together suggests that highly connected structures may be affected to the same degree by HIV infection.

All these findings require further investigation to establish whether they represent true effects of PHIV or are anomalous findings specific to our small sample. In addition to small sample size, other limitations of our study include the unequal size of the control and CPHIV groups and the cross‐sectional nature of the study. The ability to predict CPHIV based on neuroimaging is not clinically advantageous, particularly in resource‐limited settings, however identification of a neuroimaging signature of CPHIV may offer a complementary dimension to the clinical assessments for evaluating neurological effects in CPHIV. Future work will examine the link between neuroimaging features, including cortical volumes and resting state fMRI, and cognitive performance in these children.

CONCLUSION

The multimodal signature of cART‐treated perinatally acquired HIV in 7‐year old children is smaller right and left globus pallidus, bilateral amygdala, right vessel of the basal ganglia, optic chiasm, brainstem, right lateral ventricle, left ventral diencephalon, and right cerebellum white matter; enlargement of the ventricular system; increased RD and MD in bilateral uncinate fasciculus and increased AD in right uncinate fasciculus, MD and RD in right corticospinal tract, and MD in cingulum and RD in right cingulum, and AD in right superior temporal longitudinal fasciculus. Using logistic regression this multimodal signature has a classification performance of 83% sensitivity, 80% specificity, 81% accuracy, and 78% AUC, which is better than using volumes (71% AUC, 74% accuracy) or DTI measures (62% AUC, 72% accuracy) separately. Measures of cortical thickness, gyrification, and spectroscopy in basal ganglia, peritrigonal white matter, or midfrontal gray matter, are not useful in the classification of cART‐treated pediatric HIV.

Through cross‐validation and repeated resampling, EN is optimized to perform well on data that was never seen during training. The regularization and embedded feature selection ensure that both model complexity and information loss is minimized. As such, detected differences are more likely to be generalizable (and not driven by a subgroup) not only within the study population in general. Future work should examine how well these models perform on other CPHIV.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

ETHICS STATEMENT

Study approved by the Human Research Ethics committees of the Universities of Cape Town and Stellenbosch (HREC REF 123/2019 linked to 448/2011). Parents or guardians provided written informed consent and oral assent was provided by the children.

FIGURE S1 Thirty‐seven regional volumes from FreeSurfer segmentation that were included as candidate features in the sMRI volumes feature set. The 38th feature is CSF, the sum of lateral, 3rd, and 4th ventricle volumes. Image: ggseg package in R (Mowinckel & Vidal‐Piñeiro, 2020).

Click here for additional data file.

FIGURE S2 Thirty‐four FreeSurfer cortical regions in each hemisphere (Desikan et al., 2006) that were included in the sMRI cortical thickness and sMRI LGI candidate feature sets. Image: ggseg package in R (Mowinckel & Vidal‐Piñeiro, 2020).

Click here for additional data file.

FIGURE S3 Twenty white matter tracts defined by the JHU atlas (Mori et al, 2005). Candidate features for DTI in this study comprised FA, AD, RD, and MD in these tracts.

Click here for additional data file.

TABLE S1 List of candidate features for 1H‐MRS included in this study. Neurometabolite estimates from LCModel (Provencher, 2001) were obtained in the basal ganglia, peritrigonal white matter, and midfrontal gray matter.

TABLE S2 Relevant features of DTI ranked by the absolute value of their average weighting across 10 folds

TABLE S3 All classification performance metrics in each concatenation step.

TABLE S4 Classification performance metrics of the best performing feature sets when only the children common to all the feature sets are used in the analyses (n = 102).

TABLE S5 Selected features of the sMRI volumes + DTI feature set ranked by the absolute value of their weighting.

Click here for additional data file.