A Multivariate Pattern Analysis of Metabolic Profile in Neurologically Impaired Children and Adolescents

Background: The prevalence of pediatric metabolic syndrome is usually closely linked to overweight and obesity; however, this condition has also been described in children with disabilities. We performed a multivariate pattern analysis of metabolic profiles in neurologically impaired children and adolescents in order to reveal patterns and crucial biomarkers among highly interrelated variables. Patients and methods: We retrospectively reviewed 44 cases of patients (25M/19F, mean age 12.9 ± 8.0) with severe disabilities. Clinical and anthropometric parameters, body composition, blood pressure, and metabolic and endocrinological assessment (fasting blood glucose, insulin, total cholesterol, high-density lipoprotein cholesterol, triglycerides, glutamic oxaloacetic transaminase, glutamate pyruvate transaminase, gamma-glutamyl transpeptidase) were recorded in all patients. As a control group, we evaluated 120 healthy children and adolescents (61M/59F, mean age 12.9 ± 2.7). Results: In the univariate analysis, the children-with-disabilities group showed a more dispersed distribution, thus with higher variability of the features related to glucose metabolism and insulin resistance (IR) compared to the healthy controls. The principal component (PC1), which emerged from the PC analysis conducted on the merged dataset and characterized by these variables, was crucial in describing the differences between the children-with-disabilities group and controls. Conclusion: Children and adolescents with disabilities displayed a different metabolic profile compared to controls. Metabolic syndrome (MetS), particularly glucose metabolism and IR, is a crucial point to consider in the treatment and care of this fragile pediatric population. Early detection of the interrelated variables and intervention on these modifiable risk factors for metabolic disturbances play a central role in pediatric health and life expectancy in patients with a severe disability.

Considering the improvement in survival rates of patients with disabilities over the past half a century [22][23][24], the care and management of this group are of particular importance in order to improve comorbid metabolic issues, to promote their chances of growth and survival, and to prevent premature mortality.
Multivariate pattern analysis is a widely applied biomedical research tool used in the treatment and diagnosis of diseases [25,26] for its ability to take into account the structure of correlation and explore relevant associations within the data.
The aim of this study was to perform a multivariate pattern analysis of the metabolic profile in neurologically impaired children and adolescents in order to reveal patterns and crucial biomarkers among highly interrelated variables. Early detection of metabolic disorders may be crucial for the treatment and care of patients with a severe disability.

Patients
We retrospectively considered 44 patients (25 males, 19 females, mean age 12.9 ± 8.0 yrs standard deviation (SD)) with severe disabilities (Level 5 according to the Gross Motor Function Classification System) [27] who were bedridden and lived at home or in sheltered communities. In all patients, body mass index (BMI) < 2SD was calculated as previously described [15,16]. The patients had been referred to our Pediatric Endocrinological Unit for auxological evaluation or to the Pediatric Surgical Unit for treatment and/or management of nutrition support. Diagnoses included cerebral palsy due to hypoxic-ischemic damage (38.6%), dysmorphic syndrome (34.1%), and epileptic encephalopathy (27.3%). In 42/44 of the patients (95.4%), anticonvulsive drugs were administered (at least two of the following: phenobarbital, phenytoin, valproic acid, topiramate, lamotrigine, carbamazepine, and clonazepam). Bolus (72.7%) or continuous (27.3%) enteral feeding was performed in all subjects.
Clinical and anthropometric parameters, body composition, and biochemical and endocrinological assessments were recorded in all patients.
As a control group, we evaluated 120 healthy Caucasian children and adolescents comparable for age and sex (61 males, 59 females, mean age 12.9 ± 2.7 years SD), who were enrolled in our center as controls for other metabolic studies. All the parents or guardians provided their consent to retrospectively enroll the subjects in other studies for clinical research purposes, epidemiology, study of pathologies, and training aimed at improving knowledge, care, and prevention. The study was performed according to the Declaration of Helsinki and with the approval of the Institutional Review Board.

Anthropometric Parameters and Body Composition
Physical examination of the participants included evaluation of weight, height and body segment, lengths according to Stevenson's method, waist circumference, BMI, pubertal stage according to Marshall and Tanner [28,29], and blood pressure (BP) measurements, as previously detailed [15,16].

Biochemical and Endocrinological Parameters
Blood samples were drawn in the morning after overnight fasting. Metabolic and hormonal blood assays included fasting blood glucose (FBG), insulin, total cholesterol, highdensity lipoprotein (HDL) cholesterol, tryglicerides (TG), glutamic oxaloacetic transaminase (GOT), glutamate pyruvate transaminase (GPT), and gamma-glutamyl transpeptidase (GGT). As previously described [15,16], plasma glucose was measured using the hexokinase-G-6-phosphate dehydrogenase method (Siemens Healthcare Diagnostics, Camberley, UK) with a chemistry analyzer (Advia XPT, Siemens). An enzymatic method (Advia XPT, Siemens Healthcare Diagnostics, U.K.) was used to determine total cholesterol. HDL cholesterol was measured by means of the selective detergent method, followed by enzymatic reactions (Siemens Healthcare Diagnostics). The glycerol phosphatase oxidase method was used to measure TG concentrations (Siemens Healthcare Diagnostics, U.K.). Serum insulin was measured with a solid-phase, two-site chemiluminescent immunometric assay with an immunochemistry analyzer (Immulite 2000, Siemens Healthcare Diagnostics, U.K.). AST, ALT, and GGT were measured with a chemistry analyzer (Advia XPT, Siemens Healthcare) equipped with dedicated reagents; the transaminase assay method is based on nicotinamide adenine dinucleotide ((NAD)H monitoring by ultraviolet (UV) detection without the addition of P-5 -P. The GGT assay method is based on the transfer of the gamma-glutamyl group from L-gamma-glutamyl-3-carboxy-4-nitroaniline to the glycylglycine acceptor to yield 3-carboxy-4-nitroaniline, which is measured.
Insulin resistance was determined by means of the homeostasis model assessment for insulin resistance (HOMA-IR) using the following formula: insulin resistance = (insulin × glucose)/22.5 [30].
The trygliceride glucose index (TyG index) was calculated as the Ln[fasting triglycerides (mg/dL)×fasting plasma glucose (mg/dL)/2] as a surrogate marker of IR and predictor of diabetes [31].

Statistical Analysis
In order to explore the anthropometric and metabolic characteristics of the two groups (children with disabilities and healthy controls), a univariate statistical analysis was performed. To estimate the probability density function of the variables analyzed, a nonparametric kernel density estimation method was applied, representing function by means of violin plots, including a scatterplot and a bar chart showing mean and standard deviation values, in order to obtain a comprehensive view of the distribution. The graphs were grouped by sex. This method allowed us to identify patterns in the distributions of the variables. An examination of the summary statistics and boxplots would not highlight these patterns sufficiently.
The analyses were conducted separately in the two groups at first and subsequently in a merged dataset for the variables that the two groups shared. To measure the linear correlation between each pair of variables, Pearson correlation coefficients were computed and shown in a matrix correlation plot.
To study the potential association between the studied parameters, a principal component analysis was performed. This method allowed us to reduce the dimensionality of the dataset by projecting each data point onto the first few principal components (up to three) while preserving as much of the data variation as possible. This method consists of rotating the axes of the multivariate space of the original variables, along orthogonal directions of maximal variance (principal components PCs), and creating a new space defined by the PCs. Each of the PCs is characterized by a percentage of explained variance of the data. If a relevant amount of variance is explained by the first PCs, the projection of the variables as vectors on the subspace defined by the new axes can be useful in exploring correlation structures present in the data. If two variables are strongly correlated, they are projected close together (the correlation is positive) or conversely with maximum alignment (the correlation is negative). Otherwise, if they do not display correlation, they tend to be projected at an angle of 90 degrees. Moreover, if the subspace defined by the two principal components describes the variables well, these are projected toward a circle of radius unity, also known as the circle of correlation; otherwise, they tend to be projected close to the origin of the axes. Only the variables common to both groups were considered in the analysis.
Finally, to assess the association between the condition of disability and the anthropometric and metabolic variables accounting for their associations, a proportional odds multivariable logistic regression model using maximum likelihood estimation was employed. Non-linear effects were taken into account by modeling the restricted cubic splines (RCS) transformation for the variables. Restricted cubic splines (RCS) with 3 knots and 2 degrees of freedom transformation were considered. The statistical significance of the coefficients was tested by combined Wald ANOVA tests, both for linear and non-linear effects. A stepwise variable selection method considering the Akaike Information Criterion (AIC) was used to assess the chance of reducing model complexity in favor of its interpretability. The coefficients resulting from the final model indicated the relative change in the log of the odds of being a child with disabilities (1) compared to being a healthy control child (0) for a unit change of the corresponding model predictor. Graphs showing the conditional effects of the variables were generated by setting the values of the adjusting variables at their medians. Furthermore, in order to test whether the regression parameters were affected by the presence of possible outliers, a robust model was fitted by using the glmRob function of the R package robustbase V093-6. All the statistical analyses were conducted with R software (version 4.0.0).

Univariate Analysis
In Table 1, the clinical features of the patients and controls are reported. The groups were comparable for pubertal stages. For the children-with-disabilities group, 44 observations were available, and they showed missing values in some of the variables (HDL cholesterol = 1, systolic BP = 6, diastolic BP = 6). For the healthy control group of children, 120 observations were available and also in this case they showed missing values in some of the variables (fasting insulin = 4, HOMA-IR = 4, HDL cholesterol = 10, total cholesterol = 10, systolic BP = 1, diastolic BP = 1). Table 1 summarizes the main statistics of the variables considered in the univariate analysis for both groups of patients. Major differences between the two groups are evident for fasting insulin, HOMA-IR, fasting triglycerides, and HDL cholesterol. Figure 1 reports the distribution of the variables under examination in the two groups: children with disabilities (panel A) and healthy control children (panel B). The violin plots, including the scatterplot, show the differences between the distribution of the variables in males and females. Overall, the children-with-disabilities group showed a more widely dispersed distribution of the variables related to glucose metabolism and insulin resistance compared to the healthy controls. Specifically, there were differences between the disabled children group and the healthy control children in the biochemical variables of fasting insulin, HOMA-IR, fasting triglycerides, and HDL cholesterol, but there were no major differences between the two sexes within either group.
Concerning the children-with-disabilities group, major differences in dispersion were found in fasting TG and the TyG index, with males showing a wider distribution compared to females. In the healthy control children, females displayed more clustered observations for BMI compared to males, whereas the distribution of fasting TG and the TyG index was more dispersed in females than in males.

Multivariate Analysis
The first step of the multivariate analysis consisted of evaluating the correlation among the variables by computing the Pearson correlation coefficients for every pairwise association after normalization of the variables for their standard deviation.
The coefficients were computed separately for the two groups and subsequently in the variables of the merged group. The results were represented in matrix correlation plots.
Since pairwise associations do not take into account the joint relationships with the other variables, a thorough investigation of the association structure among the parameters was performed by exploiting the principal component analysis (PCA).
In the case of the children-with-disabilities group, the first three principal components explained the 56% variance in the data (PC1 = 22 %, PC2 = 18 %, PC3 = 16%). The projections of the original variables on the three planes are shown in Figure 2. The first axis (PC1) was mainly characterized by the weight, height, and BMI variables (which are highly correlated according to the BMI formula). The second axis (PC2) was characterized by the fasting insulin, HOMA-IR, and fasting glucose variables. The third axis (PC3) was characterized by the fasting TG, TyG index, and total cholesterol variables. It seemed that variables such as systolic and diastolic BP were not well represented by the three planes defined by the first three PCs. In fact, the vectors of their projections were close to the origin of the axes. In all of the planes defined by the PCs, it appeared that the fasting insulin, fasting glucose, and HOMA-IRvariables were positively correlated and grouped. In the plane defined by PC1 and PC2, these were uncorrelated to the weight, height, and age variables. In the plane defined by PC1 and PC3, the fasting TG and TyG index variables were positively correlated with each other, although neither correlated with the age and weight variables.
In the case of healthy control children, the first three principal components explained the 66% variance in the data (PC1 = 35%, PC2 = 18%, PC3 = 13%). The first axis (PC1) was mainly characterized by the variables relating to BMI, weight, height, age, fasting insulin, HOMA-IR, fasting glucose, systolic BP, and diastolic BP, which were positively correlated and so grouped together. The fasting TG and TyG index variables (strictly correlated by the TyG index formula) characterized the second axis (PC2), whereas the third axis (PC3) was mainly characterized by the total cholesterol and HDL cholesterol variables.
Finally, the multivariate PCA analysis on the merged group showed that the first three principal components explained the 57% variance in the data (PC1 = 23%, PC2 = 22%, PC3 = 12%). As expected, in this case, whereas the total cholesterol and HDL cholesterol variables represented the third axis (PC3), the height, systolic BP, diastolic BP, weight, BMI, and age variables were all grouped together, and they characterized the second axis (PC2) but not the first (PC1), contrary to what was observed for the healthy control group. Moreover, the first axis (PC1) was mainly characterized by the group of variables related to insulin resistance (fasting glucose, fasting insulin, fasting TG, TyG index, and HOMA-IR). These all positively correlated with each other and positively correlated with the PC1 (Figure 3). In the following phase, we analyzed the three planes defined by the PCs for distribution of all the observations of the merged dataset, stratified and labeled by the presence or absence of disabilities. As expected, in the plane defined by the PC1 with one of the other two dimensions, the observations relating to the healthy control children (absence of disability) were clearly seen to cluster together in the region of the plane where the PC1 is negative, whereas the observations relating to the group of children with disabilities were far more frequently found in the positive values of the PC1. If the PC1 had not been considered, as in the case of the plane defined by the PC2 and PC3, no clusters formed by the observations of the healthy control children group would have been identifiable. Indeed, the observations would have been dispersed as in the case of the children-with-disabilities group, underlining the importance of the PC1 in characterizing the difference between the two groups ( Figure 2). Overall, the PC1, which resulted from the principal component analysis conducted on the merged dataset and characterized by the variables related to glucose metabolism and insulin resistance, was crucial in describing the differences between the group of children with disabilities and the group of healthy control children.  PC3 (b), we note that the observations from the healthy control children group cluster together in a region of the plane where the PC1 is negative, whereas observations relating to the children with disabilities group are more dispersed in the direction of the projection of the variables related to insulin resistance and glucose metabolism. By excluding the PC1 and considering the PC2 and PC3 (c), the observations are found to be mixed, highlighting the fact that the PC1, characterized by the variables related to insulin resistance and glucose metabolism, is crucial in differentiating the two groups.

Logistic Regression Model
An additive linear multiple logistic regression model was fitted by considering the effects of the predictor variables: age, BMI, weight, HOMA-IR, TyG index, systolic BP, and diastolic BP. Non-linear effects were considered for HOMA-IR and weight, which showed statistical evidence of departure from the linear assumption. These were fitted in a final model in which the variables age, TyG index, and the restricted cubic spline transformation of weight and HOMA-IR were statistically significant (age, p = 0.0087; weight, p = 0.023 for the linear effects and p = 0.017 for the non-linear effects; HOMA-IR, p = 0.024 for the linear effects and p = 0.011 for the non-linear effects; TyG index, p = 0.044), whereas the systolic BP (p = 0.047) and diastolic BP (0.10) predictors were included as adjusting variables.
Concerning the linear effects of age and TyG index, the estimated coefficient yielded positive values according to a proportional increase of the log odds of being disabled with the increase of these variables. As regards the non-linear effects, Figure 4 provides additional details. Namely, the HOMA-IR showed an increasing effect on the log odds of being disabled until approximately 2, after which it leveled out to a more or less constant value, whereas weight showed a steadily decreasing effect with a reduction at higher values. The panel shows the three subspaces defined by the three PCs. When the subspace defined by the two principal components describes the variables well, these are projected toward a circle unity radius, also known as the circle of correlation. As can be seen in (a,b), the trygliceride glucose index (TyG index), homeostasis model assessment for insulin resistance (HOMA-IR), and fasting insulin variables are all well described by the subspace defined by the PC1 and PC2; they all correlate with each other and also positively correlate with PC1. In (c), we can note that the subspace defined by the PC2 and PC3 describes well only the total cholesterol variable. Shadowed areas represent 95% confidence intervals. The curve is derived from a multiple logistic model adjusted for all the covariates. On the Y axis, the conditional effect in terms of the log odds of being a child with disabilities is displayed, whereas, on the X axis, the values of the specific variables are shown. As can be seen in panel (a), the HOMA-IR has a major effect in discriminating between the group of children with disabilities and healthy control children at values below 2. In fact, at values greater than 2, the curve almost reaches a plateau. In panel (b), the non-linear conditional effect of the weight variable can be observed with a decrease in log odds of being a child with disabilities as weight increases, with a reduction at higher weights.
The robustness of the model was tested by fitting a model containing the predictors considered previously. In the robust model, the effects and statistical significance of all the predictors were confirmed without any major differences. Even the fast-backward selection AIC method did not reject any of the variables considered in the above model. Therefore, as expected, the multiple logistic regression model showed a significant association between the condition of disability and the factors related to glucose metabolism and insulin resistance, in accordance with the results of the PCA previously presented.

Discussion
We described a multivariate pattern analysis of metabolic profiles in neurologically impaired children and adolescents. We showed that the variables related to glucose metabolism and insulin resistance were crucial in describing the differences between the group of children with disabilities and the group of healthy control children.
MetS is a constellation of interconnected risk factors of metabolic origin, leading to CVD and diabetes and increased risk of mortality. Due to varying definitions, the prevalence of MetS in the pediatric age remains unclear [33]. In the literature, a median prevalence of 3.3% (range, 0-19.2) in the general population is reported and is significantly higher in overweight (11.9%) and obese (29.2%) subjects [13,14,33].
Although the pathogenetic mechanism of MetS is not fully understood [34][35][36][37], the interaction between obesity, IR, and inflammation is known to play a key role in its development. It is believed that IR initiates different pathogenic pathways, which increase metabolic risk and result in the full expression of MetS [35,38]. IR leads to a decrease in its effect on the suppression of hepatic glucose production in the liver [33]; additionally, hyperinsulinemia leads to increased TG production through an increase in transcription of hepatic lipogenic genes. The increase in free fatty acid delivery results in hepatic insensitivity to the inhibitory effects of insulin on very-low-density lipoprotein (VLDL) secretion and overproduction of triglyceride-rich VLDL particles [33]. Elevated BP in MetS is considered to be related to IR via mechanisms such as sympathetic nervous system activity, renal sodium retention, and smooth muscle growth [33]. IR also influences endothelial dysfunction and altered vasodilatory response [33].
In this report, we re-examined the metabolic profile from different points of view. Applying univariate tests to individual response variables has the advantage of simplicity of interpretation. On the other hand, it fails to account for the correlation in the data. In contrast, multivariate statistical techniques might more adequately capture the multidimensional metabolic pattern of health-related conditions. In our study, PCA results and biplots of Figure 2 provide a direct view of the association among the relevant variables, which can be grouped accordingly. These kinds of results are not available from basic univariable analyses and tabulations. Multivariable analysis techniques should be used extensively and consistently in pediatrics to allow results to be compared in a more robust way and to exploit the study of multivariate association structures.
We noted that, compared to the healthy controls, the children-with-disabilities group showed a more dispersed distribution of the variables related to glucose metabolism and IR, without significant differences in males and females. Multivariate analysis helps to quickly identify interconnected patterns in unfavorable metabolic profiles for pediatric populations with disabilities; the parameters related to glucose metabolism play a critical role in the difference between the groups. Even though we considered two groups comparable for age, sex, and pubertal status, the influence of pubertal timing on the impact of IR could not be excluded in children with disabilities [39]. The results confirm that IR is a precocious metabolic alteration that could influence the severity of metabolic involvement. The contribution of IR to the progression of skeletal muscle wasting and lipolysis should also be considered [40]. Screening for the presence of IR should be recommended in subjects with disabilities to prevent the development of complete MetS and the impairment of body composition.
We recognize that our study has some limitations. The number of participants is insufficient for the development of predictive models; in the future, an increased sample size is expected to extend and validate these first results on statistical and clinical grounds. We considered non-obese children with disabilities; the effects on IR, considering HOMA-IR as well as the TyG index, may suggest that muscle or central insulin sensitivity could also be taken into account; further studies on the pathogenic mechanism of IR may be useful in defining the relationship between metabolic parameters. Nutritional status and additional factors, such as familiarity for dysmetabolic diseases, and level of physical activity, could additionally be included in the analysis to better evaluate the causal links between variables.

Conclusions
A different metabolic profile is found in children and adolescents with disabilities compared to controls. MetS is a crucial point to consider in the treatment and care of this fragile pediatric population. Careful evaluation, periodic monitoring, and measuring of the metabolic state are recommended in these children in order to better define their needs and decrease the risk of morbidity and mortality. Early detection of the interrelated variables and intervention on these modifiable risk factors for metabolic disturbances play a central role in pediatric health and life expectancy in patients with a severe disability.