Multivariate approach for longitudinal analysis of brain metabolite levels from ages 5-11 years in children with perinatal HIV infection.

Treatment guidelines recommend that children with perinatal HIV infection (PHIV) initiate antiretroviral therapy (ART) early in life and remain on it lifelong. As part of a longitudinal study examining the long-term consequences of PHIV and early ART on the developing brain, 89 PHIV children and a control group of 85 HIV uninfected children (HIV-) received neuroimaging at ages 5, 7, 9 and 11 years, including single voxel magnetic resonance spectroscopy (MRS) in three brain regions, namely the basal ganglia (BG), midfrontal gray matter (MFGM) and peritrigonal white matter (PWM). We analysed age-related changes in absolute metabolite concentrations using a multivariate approach traditionally applied to ecological data, the Correlated Response Model (CRM) and compared these to results obtained from a multilevel mixed effect modelling (MMEM) approach. Both approaches produce similar outcomes in relation to HIV status and age effects on longitudinal trajectories. Both methods found similar age-related increases in both PHIV and HIV- children in almost all metabolites across regions. We found significantly elevated GPC+PCh across regions (95% CI=[0.033; 0.105] in BG; 95% CI=[0.021; 0.099] in PWM; 95% CI=[0.059; 0.137] in MFGM) and elevated mI in MFGM (95% CI=[0.131; 0.407]) among children living with PHIV compared to HIV- children; additionally the CRM model also indicated elevated mI in BG (95% CI=[0.008; 0.248]). These findings suggest persistent inflammation across the brain in young children living with HIV despite early ART initiation.

Introduction

Treatment guidelines recommend that children with perinatal HIV infection (PHIV) initiate antiretroviral therapy (ART) early in life and remain on it lifelong (Violari et al., 2008; “WHO, 7.1.4 When to start ART in children,” 2013). Despite early ART and viral suppression from a young age, children with PHIV continue to demonstrate neurocognitive deficits. Cross-sectional imaging studies report brain alterations in PHIV children throughout childhood even with treatment (Hoare et al., 2014; Musielak and Fine, 2016; Van den Hof et al., 2019). Longitudinal studies are needed to better understand whether the observed HIV-related differences persist, resolve or worsen with age.

In an ongoing neurodevelopmental study, we are repeatedly imaging a cohort of children with PHIV from the Children with HIV Early antiRetroviral (CHER) trial (Cotton et al., 2013; Violari et al., 2008) who live in Cape Town, South Africa and age-matched uninfected (HIV−) children from the same community. The goal of our longitudinal follow-on study is to examine the long-term consequences of PHIV on the developing brain in early-treated children stable on ART and virally suppressed from a young age. Children received neuroimaging at ages 5, 7, 9 and 11 years, including single voxel proton magnetic resonance spectroscopy (1H-MRS) in three brain regions, namely the basal ganglia (BG), midfrontal gray matter (MFGM) and peritrigonal white matter (PWM). MRS, which noninvasively evaluates localized brain metabolism, provides insights into brain health and development and may present a sensitive marker of metabolic events underlying, and even preceding, structural changes. The three regions studied here were selected due to being implicated in HIV infection and their role in healthy brain development.

The BG plays a role in memory processes, language and signalling pathways of the frontal lobe (Booth et al., 2007). This structure is often affected early in HIV infection and without treatment, calcification of this region is common (Belman et al., 1986; Govender et al., 2011). In the ART era, even though BG abnormalities are less severe, studies continue to report HIV-related structural and metabolite concentration differences (Mbugua et al., 2016; Randall et al., 2017).

Between 5 and 11 years, both synaptogenesis and synaptic pruning occur at different rates across the brain (Johnson, 2003; Levitt, 2003). Frontal gray matter, important in both cognition and executive functioning, is the last lobe to develop (Casey et al., 2008). Imaging studies have reported alterations across gray matter regions in PHIV children on treatment (Herting et al., 2015; Nwosu et al., 2018; Toich et al., 2018).

During this age range, white matter development is related to ongoing myelination. The most substantial changes in white matter development have occurred by 10 years of age (Bashat et al., 2005; Lebel et al., 2008; Lebel and Beaulieu, 2011), even though white matter maturation continues well into adulthood (Reiss et al., 1996). Even in the treatment era, studies consistently report HIV-related abnormalities across WM regions (Ackermann et al., 2016; Jankiewicz et al., 2017).

Metabolites of interest in proton MRS include glycerophosphocholine (GPC) and phosphocholine (PCh), myo-inositol (mI), N-acetyl aspartate (NAA), total creatine (creatine plus phosphocreatine [Cr+PCr]) and glutamate (Glu). GPC+PCh and NAA are markers for membrane and neuronal integrity, respectively, and mI measures glial proliferation (Keller et al., 2006; Zhu and Barker, 2011). Cr+PCr, involved in the Krebs cycle, is a measure of cellular energy levels and Glu is an important excitatory neurotransmitter (Keller et al., 2006; Zhu and Barker, 2011).

Few cross-sectional MRS studies have examined age-dependent regional metabolic changes from birth to adulthood, establishing that the most significant changes are in infancy (Blüml et al., 2013; Hashimoto et al., 1995; Kreis et al., 1993; Pouwels et al., 1999; van der Knaap et al., 1990). Within studies encompassing ages 5–11 years, some age-related changes in NAA levels as well as NAA/GPC+PCh or NAA/Cr+PCr ratios are reported in gray and white matter (Blüml et al., 2013; Horská et al., 2002; Kadota et al., 2001; Keller et al., 2004; Pouwels et al., 1999), with one study (Blüml et al., 2013) finding accompanying Cr+PCr levels increasing in gray matter. However, of the studies including participants between ages 5 and 11 years, the bulk of the data presented were from infants and/or adults. As such, there is a need for additional MRS studies of school-aged children. In addition, since all of the studies to date represent children in the developed world, there is a need to expand the data from developing countries outside of the USA and Europe to establish global patterns of development.

Due to ongoing myelination and synaptic restructuring occurring during this period (Johnson, 2003; Levitt, 2003) one would expect concurrent metabolite changes. As noted above, within this age range MRS studies consistently report age-related changes in NAA levels only. Possible reasons for the lack of age-related metabolite changes in the few studies to date include wide age ranges, small sample sizes, cross-sectional designs, sensitivity of MRS to motion and field inhomogeneities making its acquisition challenging, and study of metabolite ratios only. As one study reported increasing Cr+PCr over this age range, age-related metabolite changes may be masked or reduced in studies that report metabolite ratios only.

Our own cross-sectional analyses of PHIV children from the CHER cohort and uninfected controls demonstrated inconsistent altered metabolism across ages and regions. In the MFGM at age 7, children with PHIV have higher GPC+PCh levels than HIV unexposed uninfected (HU) children (unpublished). In the BG at age 9 years, children with PHIV have lower NAA and Glu than HU children (Robertson et al., 2018). At 11 years we found reduced NAA and Cr+PCr in PWM and increased GPC+PCh levels in MFGM in PHIV children than their HU peers (Graham et al., 2020). Since metabolite levels change across childhood (Holmes et al., 2017) and differently in different regions, HIV-related developmental delays may affect levels differently with age, leading to conflicting results. Here we have incorporated the measurements of these metabolites across all recorded ages and regions into a longitudinal approach to examine HIV-related changes in metabolite trajectories from 5 to 11 years.

Although longitudinal analyses provide a unique opportunity and increased power to potentially resolve some of the inconsistencies from our cross-sectional analyses, providing insights into whether observed differences represent persistent damage or developmental delay, these analyses typically present their own challenges, such as dealing with missing data at certain time points and hardware upgrades. Here we present the results from a multivariate approach traditionally applied to ecological data, the Correlated Response Model (CRM) (Hui, 2016), and compare these to results from a more traditional multilevel mixed effect modelling (MMEM) approach.

Methods

Participants

Participants were 89 PHIV children from the CHER trial and 85 age-matched uninfected controls (HIV−) who received neuroimaging at ages 5, 7, 9 and 11 years as part of our ongoing longitudinal study of brain development in children with PHIV.

On the CHER trial (Violari et al. 2008; Cotton et al., 2013), newborns with PHIV and CD4+ percentage (CD4%) ≥ 25% were randomized to receive ART before age 12 weeks with planned interruption after either 40 or 96 weeks, or to start ART on presentation of clinical symptoms of disease or when their CD4% dropped below 25% in the first year or below 20% thereafter. The first line ART regimen consisted of Zidovudine (ZDV) + Lamivudine (3TC) + Lopinavir-Ritonavir (LPV/r, Kaletra). Children initiated ART on average at age 14.5 (±13.11) weeks and all were on ART by the age of 75.7 weeks. The HIV− children were recruited from a linked vaccine trial (Madhi et al., 2010) and include uninfected children born to both HIV infected and uninfected mothers.

Children were scanned using a single-channel head coil on a 3 T Allegra (Siemens, Erlangen, Germany) at ages 5 and 7 years, and following a scanner upgrade, on a 3 T Skyra (Siemens, Erlangen, Germany) with a 32-channel head coil at ages 9 and 11 years. A subset of children (7 PHIV, 17 HIV−) who had been scanned at age 9 years on the 3T Allegra prior to the upgrade were re-scanned on the 3T Skyra following the upgrade. Not all subjects contributed data at every age; 141 children provided usable data at 2 ages, 59 at three ages, and 24 at all ages. Due to a technical issue with data acquisition at age 5 years (not related to the MRS data), a subset of children (11 PHIV, 0 HIV−) were re-scanned within 6 months, providing two MRS scans at this age. Ethical approval for this study was provided by the Human Research Ethics Committees of the Faculties of Health Sciences of both the Universities of Cape Town and Stellenbosch. Parents provided written informed consent and children assent. Fig. 1 summarizes the number of children providing some usable data at each time point, the brain regions where MRS was collected and the metabolite levels examined. Although the same cohort was being studied across ages, more participants provided usable data at older ages due to being able to remain motionless for longer, increasing familiarity with the MRI scanning procedure, and the newer 3T Skyra MRI producing more high- quality data.

Fig. 1.

Spectroscopy data were collected in three brain regions – midfrontal gray matter (MFGM), basal ganglia (BG) and peritrigonal white matter (PWM) – at 4 different ages on either 3 T Allegra or 3 T Skyra MRI scanners. Subject numbers indicate the number of participants that provided some usable data at each time point. A subset of participants (n=11 PHIV; n=0 HIV−) were scanned twice (on average 7.6 months apart) at age 5 years – both times on the Allegra. At age 9 years, a subset (n=7 PHIV; n=17 HIV−) were scanned twice – first on the Allegra and again on the Skyra with an average time of 10.8 months between scans.

Data acquisition and processing

The scanning protocol included a high-resolution T1-weighted 3D EPI-navigated (Tisdall et al., 2012) multi-echo magnetization-prepared rapid gradient echo (MEMPRAGE) (van der Kouwe et al., 2008) acquisition for voxel placement, followed by three single voxel 1H-MRS point resolved spectroscopy (PRESS; 1.5 × 1.5 × 1.5 cm3 voxel; TR 2000 ms, TE 30 ms, 64 averages, 2:16 min, 4 preparation scans) acquisitions in the BG, MFGM and PWM, respectively. At the younger ages on the Allegra, an echoplanar imaging (EPI) volumetric navigated (vNav) PRESS (Hess et al., 2011) sequence with prospective motion- and shim correction was used to improve spectral quality (Hess et al., 2013), bandwidth was 1000 Hz and vector size 512. On the Skyra, bandwidth was 1300 Hz and vector size 1024. A single water reference scan (TR 4000 ms, TE 30 ms, 2 preparation scans) was obtained in each region without water suppression. Voxel positioning was performed visually by matching anatomical landmarks as shown in Figure 1.

Individual spectra from the Allegra were frequency and phase-corrected offline before averaging (Hess et al., 2014). Skyra data was averaged on the scanner. LCModel (version 6.3-1) (Provencher, 2001) was used for eddy current correction using the water reference data (Klose, 1990) and spectral quantification. LCModel zero-fills the time-domain data by a factor of two, and automatically performs phasing, referencing and quantitation on the frequency domain spectra. The LCModel basis set used for model fitting included L-alanine, aspartate, creatine (Cr), negative creatine methylene (−CrCH2), g-aminobutyric acid (GABA), glucose, glycerophosphocholine (GPC), glutamine (Gln), glutamate (Glu), myo-inositol (mI), lactate, N-acetylaspartate (NAA), N-acetylaspartylglutamate (NAAG), phosphocholine (PCh), phosphocreatine (PCr), scyllo-inositol, and taurine. Macromolecule and lipid resonances were modelled using additional Gaussian basis functions (5 macromolecules, 4 lipids) (LCModel User’s Manual, n.d.). Example spectra can be found in the appendix in Figure a.3.

SPM12 (http://www.fil.ion.ucl.ac.uk/spm) was used to segment the voxel into gray matter (GM), white matter (WM) and cerebrospinal fluid (CSF) for partial volume correction. The water concentration in the voxel was calculated using an assumed water concentration of 55.5 mol/l for CSF, 43.3 mol/l for GM, and 35.9 mol/l for WM. The voxel water concentration based on the fractional tissue composition was input to LCModel to obtain absolute metabolite concentrations of GPC+PCh, mI, NAA, Cr+PCr and Glu by scaling to the measured water reference signal (Gasparovic et al., 2006). Metabolites with percent standard deviation greater than 25% were excluded from all analyses.

The amount of GM and WM in the spectroscopy voxel was controlled for by including the relative as a covariate. Children’s ages were re-scaled by subtracting 5, so that age effects were measured from the start of the observation period and not from birth. We also controlled for different scanners.

Longitudinal Analyses

Since our data contains repeated measures for five brain metabolites across time (within subjects) and within different regions at each time point we have multiple correlated responses for each subject. We present the results based on a multivariate approach traditionally applied to ecological data, the Correlated Response Model (CRM) (Hui, 2016). A comparison to the more traditional MMEM is included in the supplementary material.

CRMs combine multivariate generalized linear mixed effect models (GLMM) with a latent factor analytic model. Adding the latent factor component for factor analysis (FA) on the residuals of the model (Hui, 2016) captures associations between the multiple response variables not already accounted for by the GLMM part of the model. This model fits separate GLMMs to each column (response variable) of the multivariate response matrix and allows for residual correlation between these columns through latent variables (Hui and Blanchard, 2018). The GLMM models include subject-specific random effects to account for the correlation induced by the repeated measurements over time for each subject.

We thus jointly fit a GLMM with a factor analytics (FA) component which results in a factor analysis of residuals as opposed to performing a factor analysis on the response matrix pre model fitting. We chose to include three latent variables in this latent factor model and all subsequent CRMs. Three latent variables were included for ease of plotting and interpretation and because 67% of the total variability was explained by the first (27.3%), second (20.5%) and third (19.6%) factors.

The model formulation for each of the 5 metabolites in each of the 3 regions is as follows: where: in which

Y = the concentration (in mM) for metabolite-region combination

j = 1, .., 15 (5 metabolites in each of 3 regions) at occasion

t = 1, …, n for subject i = 1, … 163,

β0 = response-specific offset to account for different scales (for the respective metabolite-region combinations),

X1 = the age of the child at each scan occasion,

X2 = the scanner the measurements were performed on (0=Allegra, 1=Skyra),

X3 = the HIV status of the child (1=PHIV, 0=HIV−),

X4, X5, X6 = the relative percentage of GM in the voxel for each measurement in the BG, PWM and MFGM, respectively,

X1, * X3 = interaction between HIV status and age,

u = the random effect term for repeated subject observations, and is assumed to follow a N(0; σ) distribution,

z (n × k) = the k latent variables included in the model, induced by the between response associations, and

θ (k × 15), = the loadings of the measured variables on the k latent constructs.

The loadings of the latent component of the CRM contain information about the association between the measured variables and the latent factors. These associations, together with association between observations (subjects) can be illustrated by creating biplots of observations and variables in the dimensions determined by the latent variables.

The CRM provides a GLMM for each response variable and is thus a true multivariate model. To fit the CRM models, we used the boral library in R (Hui, 2016; Hui and Blanchard, 2018) which uses Bayesian estimation procedures. Thus, for the CRM results, CI refers not to confidence intervals, but instead credibility intervals which we will refer to as CI*.

Model validity was assessed by looking at a normal quantile plot to assess the normality of the data. The Dunn-Smyth residuals (Dunn and Smyth, 1996) were plotted against linear predictors, column index and row index to identify any structures in the data that would question the validity of the model. Model fit was assessed by looking at the residual correlation from each model and the relative reduction in residual correlation due to model changes. A positive value indicates a decrease in residuals and hence an improved model fit. Anonymised data can be shared with other researchers once a data sharing agreement has been completed.

Results

Table 1 presents sample demographics of children that provided usable data at each time-point. Groups did not differ on age nor sex at any of the time points (all p-values > 0.27) and show similar gender and age distributions for PHIV and HIV− children at the different ages.

Table 1

Sample Characteristics.

Age Group	HIV Status	Sex (F)	Age (Mean ± Standard Deviation)	Scanner Used
Age 5 (n=83)[1]	PHIV (n=53)	31 (58%)	5.40±0.31	Allegra
	HIV− (n=30)	14 (47%)	5.60±0.43
Age 7 (n=130)	PHIV (n=73)	37 (51%)	7.21±0.13	Allegra
	HIV− (n=57)	25 (44%)	7.23±0.13
Age 9 (n=125)[2]	PHIV (n=74)	38 (51%)	9.25±0.24	Allegra + Skyra
	HIV− (n=51)	23 (45%)	9.37±0.61
Age 11 (n=137)	PHIV (n=77)	40 (52%)	11.63±0.27	Skyra
	HIV− (n=60)	25 (42%)	11.57±0.25

11 (11 PHIV, 0 HIV−) of these children were scanned twice at age 5 years.

24 (7 PHIV, 17 HIV−) of these children provided two sets of observations for the 9-year scan – one using the Allegra scanner and the other using the Skyra.

Note: not all children had successful acquisitions for all regions and ages due to movement and/or non-compliance.

Fig. 2 illustrates the observed metabolite profiles for children with PHIV compared to HIV− controls. The figure shows increasing linear trends were not influenced by HIV status. The exceptions are GPC+PCh, where the children living with PHIV had higher levels than HIV− children across regions, and NAA in the BG region where the slope for children living with PHIV was different compared to HIV− controls.

Fig. 2.

Longitudinal profiles of metabolite concentrations in basal ganglia (BG), midfrontal gray matter (MFGM) and peritrigonal white matter (PWM) regions for children with PHIV (blue) and uninfected controls (pink) between the ages of five and eleven.

The results for both the GLMM and latent factor (LF) components of the CRM model are given in Appendix A.1. Estimates for the effect of age and HIV infection status are presented in Tables 2 and 3, respectively.

Table 2

The effect of age on metabolites concentrations estimated by the CRM longitudinal analysis.

Region	Metabolite	Correlated Response Model (CRM)
		Estimate	95% CI *
BG	GPC+PCh	0.015	[0.003; 0.025]
	mI	0.133	[0.076; 0.182]
	NAA (age)[1]	0.042	[−0.014; 0.104]
	NAA (age)[2]	0.094	[0.025; 0.160]
	NAA (HIV* age)[3]	−0.077	[−0.137; −0.027]
	Cr+PCr	0.178	[0.109; 0.250]
	Glu	0.137	[0.029; 0.249]
MFGM	GPC+PCh	0.036	[0.022; 0.049]
	mI	0.110	[0.044; 0.176]
	NAA	0.075	[−0.004; 0.158]
	Cr+PCr	0.118	[0.050; 0.182]
	Glu	0.098	[−0.017; 0.241]
PWM	GPC+PCh	0.007	[−0.008; 0.021]
	mI	0.047	[−0.008; 0.107]
	NAA	0.228	[0.161; 0.299]
	Cr+PCr	0.097	[0.057; 0.132]
	Glu	0.133	[0.039; 0.217]

As Bayesian estimation processes were used these are credibility intervals instead of confidence intervals.

The estimate from the model without the HIV* age interaction term.

The estimate from the model with the HIV* age interaction term.

The estimate of the HIV* age interaction term.

Table 3

The effect of HIV status on metabolite concentrations estimated by the CRM longitudinal analysis.

Region	Metabolite	Correlated Response Model (CRM)
		Estimate	95% CI *
BG	GPC+PCh	0.067	[0.033; 0.105]
	mI	0.122	[0.008; 0.248]
	NAA (HIV)[1]	0.065	[−0.052;0.199]
	NAA (HIV)[2]	0.381	[0.131; 0.628]
	NAA (HIV* age)[3]	−0.077	[−0.137; −0.027]
	Cr+PCr	0.081	[−0.080; 0.234]
	Glu	0.059	[−0.155; 0.302]
MFGM	GPC+PCh	0.096	[0.059; 0.137]
	mI	0.261	[0.131; 0.407]
	NAA	0.085	[−0.077; 0.284]
	Cr+PCr	0.119	[−0.014; 0.271]
	Glu	−0.038	[−0.338; 0.247]
PWM	GPC+PCh	0.061	[0.021; 0.099]
	mI	0.035	[−0.110; 0.146]
	NAA	0.036	[−0.125; 0.191]
	Cr+PCr	−0.012	[−0.086; 0.069]
	Glu	−0.118	[−0.312; 0.064]

As Bayesian estimation processes were used these are credibility intervals instead of confidence intervals.

The estimate from the model without the HIV* age interaction term.

The estimate from the model with the HIV* age interaction term.

The estimate of the HIV* age interaction term.

The age*HIV status interaction term identified a significant difference in NAA levels within the BG region and not for other responses. Accordingly, the results within Tables 2 and 3 were obtained from the simpler model without the interaction term, with the estimates from the more complex model presented only for the effect on BG-NAA.

When looking at the effect of age, estimates in Table 2 show a significant positive effect of age on metabolite concentrations for all metabolites, with the exception of BG-NAA, MFGM-NAA, MFGM-Glu, PWM-(GPC+PCh) and PWM-mI levels, where the 95% credibility intervals included zero, thus indicating a lack of statistical significance at the 5% level.

However, the effect of age on BG-NAA is significant when we incorporate the age*HIV status interaction, giving a significant positive relationship between BG-NAA and age in HIV− children. The HIV*age interaction has a negative impact on this effect, thus, the increase in BG-NAA observed over time is diminished within children with PHIV.

The simpler models without age*HIV status interactions show the effect of HIV across all ages (Table 3). We observed significantly higher levels of GPC+PCh across all ages and regions, and higher mI in the BG and MFGM regions, in children with PHIV compared to HIV− controls. In the more complex model with age*HIV status interactions added, the estimate for HIV status shows the difference in concentrations between children with PHIV and HIV− controls specifically at age 5, while the age*HIV status interaction term shows the change in this effect with increasing age. As mentioned previously, this age*HIV interaction effect only affected the concentration of BG-NAA significantly. This finding indicates that significant difference in BG-NAA concentrations between children with PHIV and HIV− children at age 5 diminishes with age.

Fig. 3 illustrates a biplot for the first 2 dimensions of the latent model of the residuals after the GLMM component was added in the CRM model. It illustrates the association between patient observations based on residual region-specific metabolite responses and the associations among these residual metabolite responses on the same set of axes. Clusters on this biplot will reveal residual groupings between metabolite responses not explained by variables included in the GLMM model and whether they are related to HIV-status of subjects. In these plots the positions of the response variables are indicated using labelled arrows that point to the centre of the plotting positions, while the observations on subjects at different time points are indicated by coloured hollow dots stratified by HIV status. Biplots for pairwise combinations of all 3 latent axes, before and after fitting the GLMM component, are in Figure a.1 of the supplementary material.

Fig. 3.

The biplot comparing the first two latent factors for the correlated response model. Projected response variables are colour-coded by region – dark green for midfrontal gray matter (GM), orange for basal ganglia (BG), and lime green for peritrigonal white matter (WM). The pink arrows highlight the choline (GPC+PCh) response variables that cluster together independent of region.

In the biplot we see a clustering of response variables by region, except for GPC+PCh responses which cluster together (as indicated by the pink arrows) independent of region (Fig. 3). The biplot does not show groupings by HIV status, whether before or after the inclusion of HIV in the GLMM part of the model. (Fig. 3 and a.1).

The MMEM model formulation and results are given in Appendix A.2 in the Supplementary material. Tables a.4–a.7 in Appendix A.3 compare the results obtained from the two modelling approaches. The results from the two approaches supported each other. Both the MMEM and CRM models yielded positive slopes for metabolite measurements across age with similar effect size estimates and 95% intervals that overlap. The credibility intervals for the estimates from the CRM model were typically wider than the confidence intervals from the MMEM models. As a result, using a 5% level of significance, the MMEM approach showed a significant positive effect of age on metabolite concentrations for all but two metabolites: BG-NAA and PWM-(GPC+PCh), while the CRM approach additionally identified non-significant effects of age on MFGM-NAA, MFGM-Glu and PWM-mI levels. Similar to the CRM, the MMEM approach showed significantly higher GPC+PCh levels in children with PHIV than HIV− children across all ages and regions, but significantly higher mI in MFGM only, while the CRM approach additionally showed significantly higher mI levels in the BG.

Discussion

This paper presents the results from modelling the multivariate response vector of longitudinal metabolite concentrations in three brain regions. We observed significantly elevated GPC+PCh across regions and significantly higher BG-mI and MFGM-mI in children living with PHIV. These findings suggest persistent inflammation across the brain in young HIV positive children despite early ART initiation.

Longitudinal profiles

Typical longitudinal profiles

During the age range studied, the brain undergoes extensive maturation to support emerging cognitive and intellectual abilities, including synaptic reorganization and myelination. Corresponding changes occur in gray matter cortical thickness and basal ganglia volume, as well as white matter volumetric and microstructural alterations. However, there are few longitudinal studies of metabolic aspects of brain maturation in childhood (Holmes et al., 2017; Kadota et al., 2001).

The CRM results show that age had a consistently positive effect on most metabolites, apart from BG-NAA, MFGM-NAA, MFGM-Glu, PWM-(GPC+PCh) and PWM-mI, although the association of BG-NAA with age was significant in the HIV− group. This supports and extends our previous findings of age-related increases in absolute Cr+PCr in frontal gray matter, Cr+PCr, NAA and Glu in white matter, and NAA and Glu in BG (Holmes et al., 2017). Other studies using cross-sectional designs found age-related increases only in NAA (Blüml et al., 2013; Keller et al., 2004; Pouwels et al., 1999) or NAA/GPC+PCh (Horská et al., 2002) and Cr+PCr (Blüml et al., 2013) in GM, NAA/GPC+PCh in BG (Horská et al., 2002) and NAA or NAA/GPC+PCh ratios in WM (Blüml et al., 2013; Horská et al., 2002; Kadota et al., 2001; Keller et al., 2004; Pouwels et al., 1999).

Age-related changes in mI and GPC+PCh have not previously been reported, although an increase in mI with age has been observed in WM in children with HIV (Keller et al., 2004). In contrast to our findings, Kreis et al., (1993) report almost no change in absolute choline and mI in this age range, with a possible transient drop below adult values at around 5 years. The age-related increase in ml and GPC+PCh that we observe in GM-containing regions within the 5–11 year age range could therefore reflect the return from this transient drop towards adult metabolite levels.

However, it should be noted that the choice of referencing method, namely to Cr+PCr or water, may influence reported results. Because obtaining reliable absolute concentrations (referenced to water) is more complicated, most in-vivo MR spectroscopy studies of the brain report metabolite ratios to Cr+PCr. However, the interpretation of metabolite ratios is ambiguous when the concentration of Cr+PCr changes, either during development or as a result of disease. On the other hand, if water is used as an internal reference, assumptions about voxel water content and relaxation rates must be made. Estimations of these unknown parameters may be inaccurate, and the true values may also be affected by diseases, such as oedema, that alter tissue water content.

Tissue water content is also known to change during brain development. The infant brain has a high water content, which falls over the first months of life, in GM reaching almost its adult value by 5 years (300 weeks gestational age), but continuing to decrease gradually in WM during early childhood (Kreis et al., 1993). The long transverse relaxation times of water in the infant brain reach adult-like values by 5 years. In this study we assumed a constant adult-like water concentration and did not account for possible age-related changes in water concentration or relaxation, therefore the observed age-related increases in metabolite levels may largely represent a decrease in tissue water concentration with age. However, this is likely to most affect the WM voxel as Kreis et al., (1993) show very little change in GM water content or T2 relaxation between the ages of 5 and 10 years.

Because each referencing method has distinct advantages and disadvantages, it is considered good practice to present both absolute metabolite concentrations and ratios to Cr+PCr. We therefore also present age-related trajectories of metabolite ratios (Supplementary Fig. a.4), thus also enabling comparison with previous studies that present only changes in metabolite ratios in children. Unlike the absolute concentrations, within our cohort almost all metabolite ratios remain constant with age. Only in PWM do we find that NAA/Cr+PCr significantly increases with age. Previous literature consistently reports NAA increases, suggesting this is a robust finding regardless of study size or referencing method (absolute or ratio). The differing trajectories of absolute concentrations and Cr+PCr ratios means that interpretation of our results is not straightforward. The age-related increases in absolute Cr+PCr levels could point to the suppression of age-related increases in metabolite ratios or alternatively the age-related increases in absolute metabolite levels may be a consequence of alterations in tissue water with age.

HIV-related effects on longitudinal profiles

The CRM showed significantly elevated GPC+PCh levels across regions as well as MFGM-mI and BG-ml in children living with PHIV. GPC+PCh is a marker of lipid membrane breakdown and, together with mI, is often associated with inflammation and astrocytosis. In adults, elevated GPC+PCh and mI is considered an early marker of inflammation, preceding neurological injury which is reflected by NAA decrease or MRI-visible damage (Chang et al., 2003; Prado et al., 2011; Tracey et al., 1996; Young et al., 2014; Zahr et al., 2014).

Higher GPC+PCh/Cr+PCr ratios have been reported in both WM and GM in PHIV children compared to HIV− controls (Prado et al., 2011; Van Dalen et al., 2016), in GM in those with previous Centers for Disease Control and Prevention (CDC) stage C disease (Van Dalen et al., 2016), and higher GPC+PCh in MFGM in those with high viral load compared to those with low HIV viral load (Keller et al., 2004). Elevated GPC+PCh in these cases indicates glial proliferation, possibly due to activation in response to ongoing neuroinflammation. Two-dimensional spectroscopy techniques have also shown that children and youth with PHIV have higher GM mI/Cr+PCr (Banakar et al., 2008). The elevated levels of GPC+PCh along with elevated mI therefore suggest inflammation despite early treatment in these children. This is consistent with observations made by Graham et al. (2020) where elevated MFGM-GPC+PCh was identified in children with PHIV at age 11. The Cr+PCr ratio data also confirms the finding of elevated GPC+PCh in PHIV, as well as to a lesser extent elevated mI in MFGM, indicating that these results are not an artefact resulting from referencing to water.

Previously, in the same cohort at 5 years of age, we found higher absolute GPC+PCh levels in the BG in children with PHIV who started ART before 12 weeks than those who initiated treatment later and the uninfected control group (Mbugua et al., 2016). Despite not seeing elevated BG-GPC+PCh cross-sectionally at 7, 9 or 11 years (Robertson et al., 2018; Graham et al., 2020), the additional power afforded by analysing the longitudinal profiles confirmed elevated GPC+PCh levels in all three regions of the brain.

We also confirmed the elevated BG-NAA levels noted at 5 years, where children with PHIV who started ART before 12 weeks had higher NAA levels than controls (Mbugua et al., 2016). Elevated BG-NAA is unexpected, as NAA levels are typically interpreted as reflecting neuronal integrity. However, the age*HIV status interaction significantly affected BG-NAA, showing that as time progressed the difference in NAA concentrations between PHIV and HIV− children decreased. This flatter trajectory, with comparable levels to HIV− children by age 11 years, raises concerns that NAA levels in PHIV children may be lower than in their uninfected peers by the time they reach adolescence.

Region effects

An interesting observation from the CRM biplot, is the clustering of response variables by region, except for GPC+PCh levels in all regions which cluster together. This result agrees with a factor analysis of cross-sectional data at age 11 that identified a GPC+PCh latent factor based on the metabolite responses (Graham et al., 2020). This suggests that apart from GPC+PCh, there is generally a stronger within-region association between different metabolites than between any metabolite levels across different regions. Correspondingly, in a subset of data presented here, we previously found strong couplings between different metabolites within the PWM and MFGM, although not in the BG (Holmes et al., 2017).

Although within-region metabolite couplings are region- and metabolite-specific, methodological causes for these must be considered before attributing them to a physiological origin. Whereas metabolite correlations between regions are based on independent measurements from different voxels, for intra-regional correlations the metabolite measurements are obtained from the same spectral acquisition and fitting procedure, and using the same water reference. However, if within-region metabolite correlations were driven by the water reference, one would expect strong within-region correlations between all metabolites, and therefore the strong correlation in GPC+PCh between regions is still noteworthy.

Furthermore, the same caveat would apply when interpreting relationships between metabolites referenced to Cr+PCr. Despite the paucity of studies investigating relationships between metabolites within and across regions, using J-difference spectroscopy, Waddell et al., 2011 similarly found that several metabolite ratios to Cr+PCr were intra-regionally correlated within the anterior cingulate and cerebellar vermis, but that only glutamate was correlated between the two regions. Van Dalen et al. (2016), however, found stronger associations between GM and WM ratios of NAA, Glu, mI and GPC+PCh to Cr+PCr than within-region associations between these metabolites in PHIV children. This was especially marked for GPC+PCh where there was a strong relationship between GM and WM levels, but without a significant relationship in GM between GPC+PCh and other metabolites, which parallels our findings. For the other metabolites, their stronger between than within-region correlations conflict with our findings, perhaps due to differences in spectral acquisition (single voxel vs spectroscopic imaging) or referencing method (to Cr+PCr vs water) and are worth further exploration.

Scanner related effects on longitudinal profiles

In longitudinal studies there are several confounding variables requiring consideration during modelling. Change of scanner hardware or software is often unavoidable, due to failure or upgrades, or the need to scan at multiple sites. We found a significant effect of scanner across all metabolites and regions (Tables a.2 and a.7), but with varying effect sizes. Thus, although metabolites measured using the Skyra scanner had higher concentrations than the Allegra, we did not identify a systematic scanner effect. The scanner and age effects were confounded since for the most part nearly all younger children (5 & 7 yrs) were scanned on the Allegra, and the older children (9 & 11 years) on the Skyra. Repeat scans using different scanners at age 9 mitigated this confounding to some extent, and the inclusion of the scanner effect in the models allowed for the partial adjustment of possible confounding of age effects due to the different scanners that successive measurements were taken on. This emphasizes the importance of accounting for different scanners in longitudinal spectroscopy analyses.

Another consideration in longitudinal studies is that using a constant voxel size may result in different voxel tissue composition over repeated scans, as brain structure and size changes during development, particularly in children. In addition, voxel placement may not be exactly repeatable from scan to scan. We therefore included for each measurement the relative percentage of gray matter in the voxel in our models. As expected, this was significant for many of the metabolite responses, due to the known differences in function and metabolic activity between gray and white matter.

Despite finding elevated GPC+PCh across regions and regionally elevated mI in children with PHIV, the trajectories of most metabolites were not different from those of HIV negative controls. This suggests that apart from NAA in the BG, the increase in metabolite levels with age follows a typical developmental trajectory for PHIV children who received early ART (< 2 years) and have been virally suppressed from a young age. However, the low variability of GPC+PCh relative to the other metabolites, and hence greater power for comparisons involving GPC+PCh, may account for the identification of HIV status effects on GPC+PCh but not metabolites that have larger variability. Future work could investigate whether GPC+PCh and mI are related to treatment and disease characteristics of the PHIV children in our cohort.

Statistical approaches

The traditional MMEM method treats each subject-occasion combination for each region-metabolite combination as separate observations, thus relying on a multilevel hierarchical model structure as illustrated in Appendix A.2. An advantage of this approach is the allowance for unbalanced designs. Additionally, the inclusion of the interaction terms allows for formal inferences with respect to the different effects of HIV status and different age profiles across metabolites and across regions. A disadvantage in this approach is that the inclusion of many higher order interactions complicates the interpretation. The analysis of the variance components of the MMEM showed that the variabilities within regions and for metabolites were not equal, confirming the selection of a stratified variance structure (see Table a.3).

In contrast, the CRM gives separate model estimates for each metabolite and region-specific response, allowing for a straightforward interpretation of the association between the covariates of interest and the response variables. An advantage of the CRM approach is additional insight into unmeasured data characteristics not included in the GLMM component. The CRM latent factor output reveals an inter-regional association between GPC+PCh concentrations and intra-regional associations between different metabolites. This is a key advantage for the CRM approach as the MMEM method would not have identified these patterns.

A disadvantage of the CRM is that it requires complete data (i.e. individuals with data for all metabolites across the three brain regions at any specific time point, rather than across time points), and consequently some data were lost. Another limitation is that there is no way to include response specific variables in this technique. Because the same variables have to be included for each observation, variables measuring the effect of %GM in all three regions have to be included for the metabolite level response from each region. Although it does not make sense to look at the effects of, for example, %GM in the BG on metabolite levels in the MFGM, we do see that these effects are non-significant, with effect sizes close to 0.

Comparison of the effect of HIV status and age on different response combinations is through the degree of overlap or separability of credibility regions. A further extension of the CRM model would permit including variable traits (in our case region and metabolite type) thus allowing for inference focused on the interaction between the effects in the GLMM and the trait variables (Hui, 2016; Hui and Blanchard, 2018). These models typically have very low power and, in our case, did not enhance the analysis.

The MMEM approach allowed formal inference of the region-metabolite interactions while the CRM did not. However, this MMEM formulation was extremely complex while the CRM formulation was easy to interpret. Hence, depending on the analysis goal, both prove to be very useful for multivariate longitudinal analysis. Both methods show promise for extending to other types of longitudinal brain imaging data which frequently show correlation between responses.

Limitations

One limitation of this longitudinal study is the possibility of inconsistent voxel positioning across time. Voxels were placed using anatomical landmarks, however variation in positioning may have occurred due to the head growth during this age range and head orientation in the scanner. Although we accounted for the effect of different scanners by including the scanner in our models it is possible that a different approach, such as scanning a phantom with known metabolite concentrations on each scanner, would have minimized the difference in metabolite level estimates between scanners. Alternatively, the use of a tool such as ComBat which identifies and removes scanner-specific variation may have also minimised the difference between estimates (Fortin et al., 2017; Johnson et al., 2007). Finally, since we did not attempt to correct our absolute metabolite estimates for possible changes in brain water content and relaxation times with age, the longitudinal profiles should be interpreted with caution. Future work should measure age-related changes in water content and relaxation in a large, longitudinal dataset in order to facilitate accurate water referencing in studies of neurometabolic brain development in childhood. Our primary result of consistently elevated GPC+PCh in PHIV is robust in that it is apparent in both absolute and ratio data, nevertheless it would also be of use to investigate possible effects of HIV infection on brain regional water content and relaxation rates. Further investigation of intra-regional metabolite correlations should also be conducted, with independent spectra and water references acquired from the same region and exploring whether relationships remain present when referencing to Cr+PCr as well as to water.