Exploring structural connectomes in children with unilateral cerebral palsy using graph theory

Abstract We explored structural brain connectomes in children with spastic unilateral cerebral palsy (uCP) and its relation to sensory‐motor function using graph theory. In 46 children with uCP (mean age = 10 years 7 months ± 2 years 9 months; Manual Ability Classification System I = 15, II = 16, III = 15) we assessed upper limb somatosensory and motor function. We collected multi‐shell diffusion‐weighted, T1‐weighted and T2‐FLAIR MRI and identified the corticospinal tract (CST) wiring pattern using transcranial magnetic stimulation. Structural connectomes were constructed using Virtual Brain Grafting‐modified FreeSurfer parcellations and multi‐shell multi‐tissue constrained spherical deconvolution‐based anatomically‐constrained tractography. Graph metrics (characteristic path length, global/local efficiency and clustering coefficient) of the whole brain, the ipsilesional/contralesional hemisphere, and the full/ipsilesional/contralesional sensory‐motor network were compared between lesion types (periventricular white matter (PWM) = 28, cortical and deep gray matter (CDGM) = 18) and CST‐wiring patterns (ipsilateral = 14, bilateral = 14, contralateral = 12, unknown = 6) using ANCOVA with age as covariate. Using elastic‐net regularized regression we investigated how graph metrics, lesion volume, lesion type, CST‐wiring pattern and age predicted sensory‐motor function. In both the whole brain and subnetworks, we observed a hyperconnectivity pattern in children with CDGM‐lesions compared with PWM‐lesions, with higher clustering coefficient (p = [<.001–.047], ηp2=[0.09–0.27]), characteristic path length (p = .003, ηp2=0.19) and local efficiency (p = [.001–.02], ηp2=[0.11–0.21]), and a lower global efficiency with age (p = [.01–.04], ηp2=[0.09–0.15]). No differences were found between CST‐wiring groups. Overall, good predictions of sensory‐motor function were obtained with elastic‐net regression (R 2 = .40–.87). CST‐wiring pattern was the strongest predictor for motor function. For somatosensory function, all independent variables contributed equally to the model. In conclusion, we demonstrated the potential of structural connectomics in understanding disease severity and brain development in children with uCP.


| INTRODUCTION
The complexity of the human brain is thought-provoking and still not fully understood. The human brain involves multiple, closely intertwined, networks that allow us to perform highly skilled movements and cognitive processes. Advancements in Diffusion magnetic resonance imaging (dMRI) and tractography techniques have enabled the investigation of both micro-and macro-structural changes in the brain's white matter, enriching our understanding of the complex brain architecture and its development.
A brain lesion occurring early in life can disrupt the typical development of these brain networks and even result in alteration of the brain's anatomical architecture, as is the case for children having cerebral palsy (CP). CP is the most common cause of disability in children (Oskoui et al., 2013), and is defined as a group of permanent disorders of the development of movement and posture, causing activity limitation, that are attributed to non-progressive disturbances that occurred in the developing fetal or infant brain; the motor disorders of CP are often accompanied by disturbances of sensation, perception, cognition, communication, and behaviour, by epilepsy, and by secondary musculoskeletal problems (Rosenbaum et al., 2007). In one third of these children, sensory-motor impairments are predominantly present on one side of the body, also referred to as unilateral cerebral palsy (uCP) (Surveillance of Cerebral Palsy in Europe, 2000).
Multiple studies have used Diffusion MR imaging in an attempt to improve our understanding of structure-function relationships in children with uCP, in particular regarding upper limb functionality, indicating a relation between the severity of the underlying brain damage and upper limb sensory-motor function (Mailleux et al., 2021;Mailleux, Franki, et al., 2020). However, most studies focused on specific regions of interest, mostly part of the primary sensory-motor areas, while motor actions go beyond activity in the primary sensory-motor network and involve multiple networks across the brain (Mailleux, Franki, et al., 2020). Moreover, compared with typically developing children, studies in CP using whole brain dMRI structural connectivity analyses revealed a reduced white matter connectivity not only in sensory-motor regions, but also in nonmotor areas (Arrigoni et al., 2016;Ballester-Plané et al., 2017). This underlines the need for further study of the brain structural networks which would allow us to investigate how the structural network reorganizes following brain damage early in life and how this relates to sensory-motor function.
Indeed, we could expect that in the case of CP, brain changes occur on a network scale through individual neuroplastic processes.
Graph theory (GT) analysis can investigate the brain network characteristics using dMRI tractography, providing information on both global and local integration in the brain (Sporns, 2013;Sporns, 2018).
In addition, it can reveal changes in brain organization that extend beyond the injured brain areas. Such information is essential for understanding global and long-term functional outcomes of focal structural brain injury. Moreover, contrary to other analysis techniques, GT-analyses of structural connectomes are not dependent on a single region or structure that may be distorted or absent due to severe brain pathology in these children. However, the construction of structural connectome does require an accurate parcellation map from which to derive the nodes of the network. Yet, brain parcellation is challenging in populations with large lesions and notable distortion of the normal brain structure, as often encountered in children with CP. Therefore, we applied a lesion inpainting method (Radwan et al., 2021) to minimize parcellation errors, which allowed us to investigate the lesioned hemisphere and include patients with large lesions who would have otherwise been excluded, without interpreter bias.
Additionally, we used a tailored processing pipeline that accounted for the presence of lesions in all analysis stages.
In this study, we will explore the structural connectomes of the dominant (contralesional) and nondominant (ipsilesional) hemisphere in children with uCP, using graph theory analysis. We will first investigate whether structural brain connectomes differ across uCP-specific classification groups (i.e., lesion type and corticospinal tract (CST) wiring pattern). Regarding lesion type, we hypothesized a more damaged structural brain network in children with late-onset lesions predominantly affecting grey matter structures (i.e. cortical and deep grey matter, CDGM) compared with children with early-onset lesions predominantly affecting the white matter (i.e., periventricular white matter, PWM). Regarding the CST-wiring pattern, we hypothesized that children in whom ipsilateral CST-projections control the impaired hand (i.e., ipsilateral and bilateral CST-wiring pattern), would have a more damaged structural brain network, compared with children with a contralateral CST-wiring pattern. Secondly, we will explore the added value of GT-measures in predicting upper limb sensory-motor function along with other neurological factors (i.e., lesion type, CSTwiring pattern and lesion volume).

| Participants
Children with a predominantly spastic type of uCP, aged 5 to 15 years were recruited between May 2014 and April 2017 via the CP-care program of the University Hospitals Leuven (Belgium). Exclusion criteria were: (1) botulinum toxin-A injections 6 months prior to testing, The study was approved by the Ethical Committee of the University Hospitals Leuven (S55555, S56513) and parental written informed consent was obtained for all children prior to participation, according to the Declaration of Helsinki. Children aged 12 years or above were additionally asked for their written assent prior to participation.

| Clinical assessment
The evaluation of sensory-motor impairments included grip force and somatosensory function (i.e., two-point discrimination and stereognosis). Grip force was evaluated with the Jamar dynamometer (Lafayette Instrument Company, Lafayette, IN, USA), using the mean of three maximum contractions of each hand. Two-point discrimination and stereognosis were assessed according to Klingels et al. (2010) Briefly, two-point discrimination was examined distally at the index finger using an aesthesiometer to identify the minimal distance at which one or two points could be correctly distinguished. Stereognosis was evaluated via tactile identification of six familiar objects.
At activity level, bimanual performance was assessed using the Assisting Hand Assessment (AHA) (Holmefur & Krumlinde-Sundholm, 2016;Krumlinde-Sundholm et al., 2007). During a videorecorded semi-structured play session, the AHA evaluates the spontaneous use of the impaired hand during bimanual activities. Afterward, 22 items were scored and converted to 0-100 logit-based AHA units. Unimanual capacity was assessed at both hands with the Jebsen-Taylor Hand Function Test (JTHFT), evaluating movement duration during the execution of six unimanual tasks (Araneda et al., 2019).

| MRI acquisition
All MR images were acquired using the same scanner (3 T Philips Achieva, 32-channel phased-array head coil) and scanning protocol.
To limit motion, a familiarization protocol was used with children up to 10 years of age (Verly et al., 2018), and all children were able to watch a movie during scanning. Multi-shell diffusion-weighted images were acquired with spatial resolution = 2.5 Â 2.5 Â 2.

| Image processing
All images were semi-automatically processed by combining publicly available toolboxes and in-house developed scripts (using Bash and Matlab v2020b). Below we briefly discuss the processing steps and give a visual overview of the processing pipeline ( Figure 1). A detailed description is given in Data S1 methods.
Structural connectomes were constructed with MRTrix3 for each subject with TCK2connectome (Smith et al., 2015) using the whole brain, SIFT2-weighted, tractograms and modified Desikan-Killiany parcellations (Desikan et al., 2006). GT-analysis was performed using inhouse developed MATLAB scripts and the Brain Connectivity toolbox (v2019-03-03). GT-measures of characteristic path length, global and local efficiency and clustering coefficient were calculated on multiple levels, namely: the whole brain (using all 88 nodes), separate ipsilesional (using all 41 supratentorial ipsilesional nodes) and contralesional (using all 41 supratentorial contralesional nodes) hemispheres, as well as the full sensory-motor network (SMN, 22 nodes, see Table A in Data S1), and ipsilesional and contralesional SMN (each 11 nodes, see Table A in Data S1). Self-connections were removed from all connectomes. Nodes where less than 1000 streamlines arrived were defined as disconnected nodes and were removed from the connectome. Subsequent, disconnected nodes in a particular subnetwork, having no connections to other nodes within the subnetwork, were removed from that subnetwork. All edge weights were divided by the maximum edge weight. Characteristic path length and global efficiency were calculated using Dijkstra's algorithm (Dijkstra, 1959), with the connection-length matrix defined by the inverse edge weights. Clustering coefficient and local efficiency measures were calculated as recommended by Wang et al. (2017). Briefly, characteristic path length measures the average distance between any two nodes in the network. Clustering coefficient expresses the tendency of a network to be organized in densely connected groups (clusters) and is characterized by the connectivity between neighboring nodes. Global efficiency refers to the mean inverse distance between two nodes in the network. Similarly, local efficiency is a measure of the mean inverse distance between the neighbors of a node, excluding the node itself. For each subnetwork, 100 random graphs were calculated by randomly permuting the edges, while keeping the connectome symmetry, the zero-weight of the self-connections and excluding graphs with disconnected nodes. All graph measures were normalized by dividing the original graph measure by the median graph measure of the equivalent random networks.

| Transcranial magnetic stimulation
We performed single-pulse TMS to identify the underlying CST-wiring

| Statistical analysis
Descriptive statistics were collected including age, sex and side of uCP. One-way analyses of covariance (ANCOVA) were used to investigate the difference in the GT-measures between lesion types (PWM and CDGM) and between CST-wiring groups (contralateral, bilateral, ipsilateral) with age as a covariate. If a nonsignificant interaction was found between age and the clinical group, only the model with the main effects was retained. Normal distribution of the residuals was reviewed and confirmed for each fitted model (Lumley et al., 2002).
Next, we explored the value of adding GT-measures to other neurological factors (lesion volume, lesion type and CST-wiring pattern) in predicting upper limb sensorimotor function. This was done through elastic-net regularized regression which optimizes between ridge regression (L2) which shrinks coefficients, and LASSO (L1) regression which excludes predictors that do not add to the model (Zou & Hastie, 2005). The amount of ridge or LASSO regression is expressed by a variable alpha which ranges from 0 (only ridge regression) to 1 (only LASSO). The models were evaluated using R-squared 2 and root-mean-square error (RMSE). After taking the average of 100 cross-validation errors for each alpha and lambda combination, the highest alpha and lambda with a mean cross-validation error that fell within one standard error of the combination with the lowest mean cross-validation error was used. This resulted in the selection of the best, most parsimonious model. The continuous variables (GT measures, lesion volume and age) were standardized so that the estimates can be interpreted as effect sizes. Further, dummy variables were created for the categorical variables (lesion type and CST-wiring pattern). The R 2 was interpreted according to Cohen (Cohen, 1988) as weak (.02), moderate (.13) or substantial (.26). The effect sizes of the individual predictors were interpreted according to Cohen's jdj with values <.1 as tiny, values between .1 and .2 as very small, between .2 and .5 as small, between .5 and .8 as moderate, between .8 and 1.2 as large, between 1.2 and 2.0 as very large and >2.0 as huge (Sawilowsky, 2009). This analysis was performed with R (version 4.1.1).

| Participants
Fifty-five children with spastic uCP were included. Average age at time of the MRI assessment was 10 years and 7 months (SD 2 years and 9 months; age range 5 years 6 months to 15 years and 10 months). Regarding lesion type, one child was classified as having a cortical maldevelopment, six children had an acquired brain lesion and two children presented with a normal structural MRI scan. These children were excluded from further analyses as these small groups do  Table 1 displays an overview of the participant's characteristics. In addition, a detailed overview of the descriptive and clinical characteristics according to lesion types and CST-wiring pattern is provided in Data S1 (Tables B and C). As can be visually depicted, the large extent of the lesion in panel b led to multiple nodes in the connectomes to be disconnected, which appeared in 12/46 participants (median = 2 disconnected nodes, range= [1][2][3][4][5][6][7][8][9][10][11][12][13]) in this study. By removing these disconnected nodes from the connectome, these participants could be successfully included in all analyses.

| Graph theory measures across lesion type and CST-wiring pattern groups
For lesion type, significant main effects were found with higher values in the CDGM group compared with the PWM group with effect sizes ranging from medium to large (Table 2) for following parameters. Clustering coefficient was significantly higher in the whole brain (p = .048; η 2 p = 0.09); ipsilesional hemisphere (p < .001, η 2 p = 0.27), full SMN (p = .001, η 2 p = 0.21), and ipsilesional SMN (p = .01, η 2 p = 0.14). A higher characteristic path length was found in the full SMN (p = .003, η 2 p = 0.19). Local efficiency was significantly higher in the whole brain (p = .01, η 2 p = 0.15), full SMN (p = .001, η 2 p = 0.21) and contralesional SMN (p = .02, η 2 p = 0.11). Additionally, a significant interaction effect between age and lesion type was found for global efficiency in the whole brain (p = .01, η 2 p = 0.15), contralesional hemisphere (p = .03, η 2 p = 0.11) and contralesional SMN (p = .04, η 2 p = 0.09), indicating that global efficiency decreases more with age in children with CDGM lesions compared with children with PWM lesions. Another interaction effect with age was found for characteristic path length in the whole brain (p = .01, η 2 p = 0.15), which increased more with age in children with CDGM lesions. The interaction effects are visualized in Figures S1 to S5. No significant differences were found across CST-wiring pattern groups (p > .08, η 2 p <0.13), except for an interaction effect between age and CST-wiring pattern for characteristic path length in the contralesional hemisphere (p = 0.006, η 2 p =0.26) and in the ipsilesional SMN (p = .02, η 2 p =0.20), without significant post hoc differences (see Table 3 and Figures S6 and S7).

| Elastic-net regularized regression
For this analysis, 39 children with the complete data set were included. One child had missing data for two-point discrimination and stereognosis and six children did not undergo the TMS assessment.
The elastic-net regularized regression selected a full LASSOregression (α = 1) for four dependent variables (i.e., AHA, grip force of both hands and JTHFT of the dominant hand), while for both somatosensory outcomes a full RIDGE-regression was selected (α = 0). For the JTHFT of the impaired hand a combination of LASSO and RIDGE regression was used (α = .17). R 2 was substantial for all models, except for movement duration (i.e., JTHFT). A detailed overview of the output can be found in Table D in Data S1.
For bimanual performance, that is, AHA, the R 2 was .73 (RMSE = 0.56). The CST-wiring pattern was the strongest predictor with a moderate (d = À0.78) to large (d = À0.98) effect size. Compared with children with a contralateral CST-wiring pattern, having an ipsilateral CST-wiring pattern was related to a 0.98xSD lower score on the AHA on average and having a bilateral CST-wiring pattern with a 0.78xSD lower score on average. The other retained variables contributed with a tiny (jdj < 0.1) to very small (jdj = 0.1-0.2) effect.
The R 2 for grip force of the impaired hand was .74 (RMSE = 0.54).

| DISCUSSION
In this exploratory study, we used graph theory to explore structural brain connectomes across the whole brain in children with uCP, and investigated the relation with upper limb sensory-motor function. We found that structural connectomes were mainly lesion type A previous study showed that in children with uCP, CDGM lesions damages more brain regions and are more extended compared to PWM lesions (Mailleux et al., 2017). In addition, white matter tracts, such as the CST and medial lemniscus, appear to be more severely damaged in children with CDGM lesions compared with children with PWM lesions (Mailleux, Simon-Martinez, et al., 2020). Our findings additionally suggest that children with CDGM lesions have a hyperconnectivity pattern between neighboring nodes (i.e., increased clustering coefficient) in the ipsilesional hemisphere and SMN which coincided with a decreased capacity in the full SMN connectome to communicate between remote nodes (i.e., increased characteristic path length) compared with children with PWM lesions. CDGM lesions predominantly affect the target sites of the white matter connections requiring those connections to take a detour to communicate between nodes, explaining the increased characteristic path length.
Subsequently, as more direct connections between remote nodes are damaged, neighboring nodes might tend to remain connected to accommodate alternative pathways. This hyperconnectivity pattern also fits within the neural group selection theory that indicates that children with brain damage experience difficulties with selecting the most appropriate neurons to perform a motor task, and will rather select different groups of neurons every trial they perform the same motor task (Hadders-Algra, 2000). This results in the retention of abundant connections rather than pruning toward the most optimal solution. Hence, based on our findings, we hypothesize that maturational processes of pruning and myelination occur more aberrantly in children with CDGM lesions compared with PWM lesions which appears to persist into childhood and adolescence.
This hypothesis is strengthened by the interaction effects found between lesion type and age for some graph metrics. With age, global efficiency in the whole brain and in the contralesional hemisphere and SMN decreased more in children with CDGM lesions compared with PWM lesions, while characteristic path length in the whole brain increased more in children with CDGM lesions. In neurotypical development, global efficiency increases shortly after birth while both clustering coefficient and characteristic path length decreases, a process that continues even into adulthood (Cao et al., 2016). Hence, we hypothesize that the structural connectome of children with CDGM lesions deteriorates more with age compared with children with PWM lesions. Longitudinal studies will be needed to confirm this hypothesis.
So far, only one other research group has used GT-analysis to investigate structural connectomes in children with uCP following perinatal stroke, but only addressed the contralesional hemisphere (Craig et al., 2020;Craig et al., 2022). Similar to our findings, they found a higher global and local efficiency of the contralesional hemisphere and higher clustering coefficient of the contralesional SMN in children with arterial ischemic stroke (mainly affecting CDGM structures) compared with children with periventricular venous infarction (mainly affecting the white matter) as well as in both patient groups compared with neurotypical developing peers (Craig et al., 2020;Craig et al., 2022). Nevertheless, since up to 50% of children with uCP have bilateral brains lesions (Mailleux et al., 2017), future research is required to investigate to what extent these changes in the contralesional hemisphere are the result of such bilateral brain damage or rather reflect potential compensatory mechanisms.
Strikingly, only limited differences in structural connectomes were found between CST-wiring patterns. This was an unexpected finding, since reorganization of the CST is the most well-known example of brain plasticity in children with uCP. CST-projections descend from the brain to the spinal cord both contralateral as well as ipsilateral (Eyre, 2007). In neurotypical development the ipsilateral projections are gradually withdrawn, while the presence of a (unilateral) lesion may cause the contralateral CST-projections to withdraw and the ipsilateral projections to preserve (Eyre, 2007). Nevertheless, our findings suggest that the type of CST-wiring pattern occurs independently of how the rest of the brain is structurally reorganized in children with uCP, and that the structural connectome is more determined by the type of the lesion. On the other hand, in a previous publication on this cohort (Mailleux, Simon-Martinez, et al., 2020) This study also warrants some critical reflections. First, we this dataset did not include a typical developing group. However, our results, at least for the contralesional side, are in line with the studies by Craig et al. (2020); Craig et al. (2022) who additionally included a control group. Second, the heterogeneity of the included lesions and effects of motion in this pediatric cohort is challenging for conducting diffusion MRI analyses. We coped with these issues by developing a tailored processing pipeline, minimizing potential biases. Third, the acquired diffusion dataset did not include volumes with reversed phase-encoding for distortion correction. Therefore, the BDP toolbox was used to correct the EPI distortion artifacts based on the T1-weighted images. Next, the elastic-net regularized regression was performed on the whole data set, increasing the risk of overfitting the models. Nevertheless, the aim of this research question was to exploratively identify predictors. Finally, we did not correct for multiple testing due to the explorative nature of our study. The number of multiple comparisons were kept minimal by only focusing on four graph metrics across six network levels. Furthermore, effect sizes were additionally reported and strengthened the statistical findings. Moreover, the results and hypotheses made in this exploratory study could support future studies.
To the best of our knowledge, this is the first study that explored the structural connectome using GT-analysis of both the contralesional and ipsilesional hemisphere in children with uCP using a semiautomated analysis. A major strength of this study is that our diffusion MRI sequence allowed for multi-shell multi-tissue CSD-analysis which accommodates the modeling of crossing fiber populations and contributions from different brain tissues (Jeurissen et al., 2014).
Moreover, when studying the brain in children with CP using tractography, studies often lose data due to extensive lesions which then results in an underrepresentation of those children. Here, we did not lose any participant data due to the inability of tracking specific tracts owed to the use of VBG (Radwan et al., 2021) and the graph theoretical framework accounting for the effects of lesions.
In conclusion, our study demonstrated the feasibility of an automated GT-analysis in children with uCP, including children with large lesions. Our results suggest a hyperconnectivity pattern between neighboring nodes in the ipsilesional hemisphere and SMN. Additionally, the structural connectome did not differ between CST-wiring patterns. However, the CST-wiring pattern outweighed structural connectomes in predicting upper limb motor function, underlining the importance to include this variable when studying structure-function relation for the upper limb in children with uCP. For somatosensation, we could not identify a strong individual predictor. Nevertheless, graph theory analysis seems to be a powerful research tool to strengthen our insights regarding the impact of brain damage on both structural and functional connectomes of the developing brain in children with uCP and how this relates to function, as well as to capture brain changes after intensive therapy models.