Causal modelling demonstrates metabolic power is largely affected by gait kinematics and motor control in children with cerebral palsy

Metabolic power (net energy consumed while walking per unit time) is, on average, two-to-three times greater in children with cerebral palsy (CP) than their typically developing peers, contributing to greater physical fatigue, lower levels of physical activity and greater risk of cardiovascular disease. The goal of this study was to identify the causal effects of clinical factors that may contribute to high metabolic power demand in children with CP. We included children who 1) visited Gillette Children’s Specialty Healthcare for a quantitative gait assessment after the year 2000, 2) were formally diagnosed with CP, 3) were classified as level I-III under the Gross Motor Function Classification System and 4) were 18 years old or younger. We created a structural causal model that specified the assumed relationships of a child’s gait pattern (i.e., gait deviation index, GDI) and common impairments (i.e., dynamic and selective motor control, strength, and spasticity) with metabolic power. We estimated causal effects using Bayesian additive regression trees, adjusting for factors identified by the causal model. There were 2157 children who met our criteria. We found that a child’s gait pattern, as summarized by the GDI, affected metabolic power approximately twice as much as the next largest contributor. Selective motor control, dynamic motor control, and spasticity had the next largest effects. Among the factors we considered, strength had the smallest effect on metabolic power. Our results suggest that children with CP may benefit more from treatments that improve their gait pattern and motor control than treatments that improve spasticity or strength.


Introduction
The net metabolic power required during walking is, on average, two-to-three times greater in children with cerebral palsy (CP) than their typically developing peers [1,2]. As a result, (temporal priority) and all other "potential causes" (confounders) must be controlled for. Confounding is common in observational studies since data are collected in uncontrolled environments. We can reduce confounding by statistically controlling (adjusting) for factors, but we may induce bias if we fail to do so correctly. MacWilliams and colleagues have shown that without appropriate adjustment, linear regression analysis can greatly over or underestimate effect sizes for factors causing function limitations in CP [30]. If we can correctly identify confounders, causal inferences are possible [31]. Causal models provide a systematic approach to identifying confounders. One type of causal model is the structural causal model, often depicted as a directed acyclic graph. A directed acyclic graph is a graphical representation of causal relationships that can be queried to determine possible confounders [31]. Directed acyclic graphs are easier to use and less likely to fail at identifying confounders compared to traditional methods [32], enabling causal inferences from observational data.
In the study described here, we (1) propose a causal model for metabolic power in CP as a function of gait kinematics, selective motor control, dynamic motor control, strength, and spasticity, and (2) compute the total causal effect of each factor on metabolic power with large scale data. These five factors are commonly measured or treated, and likely contribute to high metabolic power in children with CP. We also included walking speed, age, height, mass, and sex within our causal model as potential confounders. We used directed acyclic graphs to form our causal model and identify confounders. Then, we used Bayesian additive regression trees (BART) to compute the causal effects of these five factors, including model-identified confounders in the analysis to adjust for their effects.

Participants
We obtained approval for this research from the Research Ethics Board (REB) at Simon Fraser University.
We retrospectively analyzed data collected from 6220 children seen between the years 2000 and 2020 at the Center for Gait and Motion Analysis at Gillette Children's Specialty Healthcare. We selected children who met the following criteria: formally diagnosed with CP; classified as level I, II, or III under the Gross Motor Function Classification System (GMFCS); 18 years or younger; and had undergone a quantitative 3D gait assessment and a 6-minute walking metabolic assessment. Many children had visited the lab more than once. To avoid pseudoreplication, we selected the visit with the least missing data.

Data collection
Spasticity was measured by a trained physical therapist using the Ashworth Scale [33,34]. The Ashworth scale is a 5-point scale defined by the following: (1) no increase in tone, (2) slight increase in tone, (3) more marked increase in tone, (4) considerable increase in tone, and (5) rigidity. Six muscle groups were assessed bilaterally: hip flexors, hip adductors, rectus femoris, hamstrings, plantarflexors, and tibialis posterior. We calculated a summary spasticity score by applying polychoric principal component analysis to the bilateral measurements from these six muscles for individuals with complete data [17,35]. We used polychoric principal component analysis since standard principal component analysis produces biased estimates with categorical data [36].
Strength was measured by a physical therapist for the hip flexors, hip adductors, rectus femoris, hamstrings, plantarflexors, and tibialis posterior. This was measured using the Kendall scale [37,38], a 5-point scale where 1 is defined as a 'visible or palpable contraction (no range of motion)' and 5 is defined as 'full range of motion against gravity'. Again, we calculated a summary strength score using polychoric principal component analysis, which included measurements from both lower limbs for individuals with complete data [17,35].
Selective motor control (SMC) was measured by a physical therapist using a 3-point scale defined by the following: (0) patterned movement, (1) partially isolated movement, and (2) completely isolated movement [14,19,35]. SMC reflects the "ability to isolate the activation of muscles in a selected pattern in response to demands of a voluntary movement or posture" [39]. The physical therapist measured SMC for hip flexion, hip adduction, knee extension, knee flexion, plantarflexion, and posterior tibialis. Again, we calculated a summary SMC score with the same methods used to calculate summary spasticity and strength scores [17,35].
Patients underwent a quantitative 3D gait analysis with electromyography of the anterior tibialis, lateral and medial hamstrings, gastroc-soleus, and rectus femoris muscles. Using the motion capture data from the gait analysis, we calculated gait deviation index (GDI) scores. We calculated a reduced order approximation of 15 kinematic gait features and scaled that value with respect to a group of children with typical development. As a result, the GDI reflects the difference in overall kinematics of an individual from the average child with typical development [40]. Using electromyographic data, we also calculated dynamic motor control (DMC) scores for each child. We calculated muscle synergies using non-negative matrix factorization and scaled the variance-accounted-for by one synergy with respect to a group of children with typical development. As a result, the DMC index reflects the complexity of motor coordination during walking [20,41].
Breath-to-breath oxygen consumption was measured during six minutes of overground, barefoot walking around a rectangular 80-meter track (Ultima CPX, Medical Graphics Corporation, St. Paul, MN, USA). Breath-to-breath oxygen consumption was also measured during 10 minutes of reclining rest preceding the walking session. The average steady state oxygen consumption was calculated for both periods according to the method proposed by [42]. Oxygen consumption per unit time (mL O 2 /s) for both periods was converted to power using the conversion rate of 20.1 Joules/mL O 2 [43,44]. We then calculated the net metabolic power by subtracting resting from walking values. We also calculated average walking speed during the trial and recorded patient age, height, mass, and sex.

Causal model
Structural causal models are visual representations of causal assumptions that we can systematically test for plausibility and query to identify confounders. Nodes represent variables, and arrows drawn (or excluded) between variables indicate causal relationships (or lack thereof) between them. The absence or presence of arrows can be used to determine conditional independencies implied by the causal model. By checking these independencies, we can test the plausibility of a model. And, just as we can test a model's plausibility, we can also determine from the model what variables are potential confounders (see Greenland & Pearl (2017) for a detailed explanation). We can then statistically adjust for these potential confounders in our computational analysis to reduce potential bias [31]. It may be necessary to adjust for more than one factor. An adjustment set is the minimum set of factors to adjust for when computing the causal effect of X on Y [31].
We define causal effects as the total effect of X on Y. This includes both the direct (unmediated) and indirect (mediated) effects of X on Y. We focused on the total effect of each factor as treatments in CP rarely affect a single factor due to their inter-relatedness. For example, strength may directly affect metabolic power, but it may also indirectly affect metabolic power via changes in kinematics.
We built a causal model for this study that represents our hypotheses and assumptions about the mechanisms contributing to metabolic demand during walking (Fig 1). We created and tested the proposed model using the dagitty package in R [45]. From our causal model, we determined adjustment sets that we used in our data analysis to compute the causal effects of each factor on metabolic power.

Data analysis
After obtaining adjustment sets from our causal model, we used the Bayesian additive regression trees (BART) method [46] to compute the causal effects of GDI, DMC, SMC, strength, and spasticity on metabolic power. Since BART is a regression-based approach, including factors from each adjustment set as predictors in the BART model adjusted for their effects and reduced confounding. A causal model identifies sources of bias, but it does not provide information regarding the strength or direction of a causal relationship. Thus, it was necessary to use regression analysis (i.e., BART) to compute effect sizes.
The BART method is an excellent option to identify the magnitude of causal effects from observational data in CP. The BART method demonstrates superior predictive abilities against other common methods [47] and has been successfully used in the causal inference of observational data [48]. Many observational studies make linear assumptions. CP is a complex condition, and although linear or multiple linear regression are convenient and easy-to-understand tools, it is unwise to assume linearity between factors. For example, as a child's tibial torsion increases (either internally or externally), we can expect their walking pattern to deteriorate, resulting in a U-shaped relationship. Besides its predictive prowess, BART is a non-parametric tool that a) uses Bayesian methods to produce estimates with honest uncertainty bounds in its predictions, b) can handle large numbers of predictor variables of various types (scale, ordinal, categorical), c) requires little-to-no hand-tuning, and d) does not require investigator involvement to determine the shape of the response surface [48]. In addition, BART natively handles missing data (increasing sample representativeness) and outperforms imputation methods when dealing with data not missing at random [49].
To visualize and interpret the causal effects specified by our BART models, we utilized accumulated local effects plots [50]. These plots show the average change in the response variable (i.e., metabolic power) due to changes in the predictor variable (i.e., GDI, DMC, SMC, strength, or spasticity). However, in conjunction with our causal model, the correct way to interpret these plots is as follows: if we were to alter one variable (e.g. DMC) by X amount while allowing all other variables affected by that variable to change naturally (e.g. GDI and speed) and adjusting for confounders (e.g. age, height, mass, etc.), we could expect metabolic power to change by Y amount. Machine learning algorithms, such as BART, are harder to interpret than linear regression because they use more complex (and often hidden) functions. By visualizing the results, we could interpret the causal effects specified by our model. Accumulated local effects plots, in particular, are unbiased when features are highly correlated compared to other visualization methods, such as partial dependence plots [50].
We used the R package bartMachine [49] to create our BART models. Following hyperparameter optimization with the function bartMachineCV, we selected the hyperparameter set with the lowest out-of-sample RMSE. As a result, our model settings were: num_trees = 50, k = 5, nu = 3, q = 0.99, use_missing_data = TRUE. For replicability, we set the seed of each model as 42. To visualize the causal effects of each factor in the form of accumulated local effects plots, we used the R package ALEPlot [51]. From each plot, we approximated effect sizes using the range of the middle 95th percentile of the sample.

Participants
We analyzed data from 2157 children who met our inclusion criteria (Table 1). Of those children, 32.1% were missing DMC scores, 18.3% were missing spasticity scores, 17.7% were missing strength scores, 17.5% were missing SMC scores, and 0.6% were missing GDI scores.

Outcomes of causal model
Our model satisfied all implied conditional independencies. The conventional cut-off for independence is ±0.3. All partial correlation coefficients were smaller than ±0.2, suggesting that the observed data was consistent with the proposed relationships in our causal model (i.e., our model was plausible). From our model, we obtained sufficient adjustment sets for all variables ( Table 2). Adjustment sets were the same for DMC, SMC, spasticity and strength.

Outcomes of BART
The purpose of this study was to determine the causal effects of GDI, DMC, SMC, strength and spasticity on metabolic power. Based on the adjustment sets indicated by our assumed causal model, we only required two BART models: one to compute the effects of GDI (r 2 = The initial brain injury is colored white to represent an unmeasured factor. Despite being unmeasured, including the initial brain injury helps to simplify the complex relationship between neurological impairments in CP. Blue nodes are factors often treated and measured in the clinic (neurological and physical impairments), so understanding their effects offers a better understanding of their importance in treating metabolic power. Metabolic power, colored in black, is the outcome of interest. Gray variables are potential confounders.
https://doi.org/10.1371/journal.pone.0285667.g001 0.81, rmse = 32.96) and one to compute the effects of DMC, SMC, strength and spasticity (r 2 = 0.72, rmse = 39.48). From the accumulated local effects plots (Fig 2), we approximated effect sizes using the range of the middle 95th percentile of the sample (Fig 3). GDI had approximately a two-fold effect on metabolic power than DMC, SMC and spasticity. Strength had the smallest effect on metabolic power. These results indicate that gait kinematics had the largest effect on high metabolic power.
The GDI had the steepest curve (Fig 2), on average, suggesting that changing GDI can elicit the largest changes in metabolic power. Most of the factors shared a near-linear or sigmoidal relationship with metabolic power. However, strength and metabolic power shared an almost inverted-U relationship. Regardless of shape, most factors plateaued at the extremes of each plot. These plateaus suggest that after a certain level of impairment, metabolic power may not increase or decrease considerably. So, children with mild or very severe level impairments may experience limited reductions in metabolic power following treatment. For example, children with largely in-tact dynamic motor control may experience limited changes in metabolic power due to ceiling effects whereas, children with severely compromised DMC may experience limited changes in metabolic power due to movement limitations.
We have reported speed, mass, age, height, and sex effects as supporting information (see S1 Appendix). Although these variables are commonly measured, they are unlikely to be the cause of high metabolic power in children with CP and thus fall outside the purpose of this study. However, we have discussed their effects below to help the reader interpret our results. Height is described in cm. Mass is described in kg. Walking speed is described in m/s. SMC, spasticity, and strength are scores derived from polychoric principal component analysis. Metabolic power is the net (walking minus resting) energy required to walk per unit time and is described in Watts. * To allow for comparison across studies, we have also provided the gross power. Gross power is the total energy required to walk per unit time, also described in Watts. https://doi.org/10.1371/journal.pone.0285667.t001

Discussion
Reducing metabolic power is a means of reducing fatigue and facilitating physical activity in children with CP. Our analysis indicates that gait kinematics (as quantified by GDI) had the largest causal effect on metabolic power, more than twice as large as the other factors. The effects of SMC, DMC and spasticity had the next largest contributions to metabolic power. Strength had the smallest effect. These results suggest that improving gait kinematics may elicit reductions in metabolic power that are more likely to be clinically meaningful. We cannot perform randomized control trials, but we must assess the order of magnitude of our results to either confirm or refute our findings. Due to the various ways that metabolic costs associated with walking are defined in the literature, assessing percent changes in metabolic power enabled the comparison of our results to experimentally determined values. On average, children with CP need to reduce their metabolic power by 50% to achieve values  within the typical range. Such a reduction may not be possible, but a 10% reduction in metabolic power is still considered clinically meaningful [52]. This value can be used to interpret the significance of the observed effects in this study. We included comparisons to distancebased measurements despite using a time-based measurement. Distance-based measures are equivalent to time-based measures when walking speed is unchanged. Since the studies discussed here either do not report, control for, or show very small changes in walking speed, we considered percent changes in distance-based measures comparable to percent changes in time-based measures.
The model-predicted changes in metabolic power following changes in GDI are comparable to experimental studies. In the study by McMulkin et al., children (GMFCS level I/II) who underwent surgical intervention with a femoral derotation osteotomy improved their GDI bỹ 13 points and their net oxygen cost (net volume of oxygen consumed per unit distance) bỹ 15% [28]. According to our results, metabolic power would decrease by 12-27W if GDI increased by 13 points. For a child with median metabolic power in our sample (124W), this would be equivalent to a 10-22% decrease in metabolic power. The values seen by McMulkin et al. fall within this range. In another study, children with CP improved their GDI by approximately five points and their net oxygen cost by 2.5% one year following single level or singleevent multilevel orthopedic surgery [29]. According to our results, metabolic power would change by 2-13% for a child with the median metabolic power if the GDI increases by five points. While our estimates seem high, Wren et al. may have found statistically non-significant results because their sample included children with more severe CP (i.e., GMFCS level IV). Children with more severe impairments (i.e., GMFCS level III and IV) have significantly greater metabolic power during walking [53]. So, even if they experience similar absolute reductions in metabolic power, relative reductions in metabolic power may not be as significant because metabolic power is higher overall. This may be why Wren et al. saw lower changes than those predicted by our results. Further, McMulkin et al. only saw significant reductions in metabolic power for children with GMFCS level I/II despite greater absolute reductions in metabolic power for children with GMFCS level III. Thus, treatments that produce greater changes in gait kinematics may yield more significant results for children with lower metabolic power when we consider relative rather than absolute changes. This research highlights how myriad factors impact metabolic power in CP, which can provide numerous pathways for potentially improving walking power. Our model suggests that changes in GDI and motor control affect metabolic power. Interventions that can improve these factors may improve metabolic power. This was recently observed in a pilot study by Conner et al. A small group of children with CP improved DMC by 7% and reduced metabolic power by 29% after resistance training with an ankle exoskeleton [23]. According to our results, metabolic power should decrease by 1-11% for a 7% improvement in DMC. These estimates are lower than Conner et al., but Conner et al. have also shown that their exoskeleton training results in mechanically more efficient gait [23] and significantly greater muscle strength [54]. Since our model does not include a metric of gait kinetics, this could have affected our causal effect estimates. In addition, a large proportion of children in our sample (31.2%) are missing DMC scores. Although BART natively handles missing data, greater missing data increases model uncertainty. Both unexplained factors and model uncertainty might explain why our estimates are lower. We estimate that SMC has similar effects on metabolic power as DMC. Clinical examination of children with CP often reveals reduced SMC due to synergistic or patterned muscle activation [39], which has been linked to inter-joint coupling during gait [55] and abnormal gait patterns [56]. We were limited to clinical examination measures of SMC, so it is possible that EMG-based measures of SMC, such as those proposed by [57], may have provided more sensitive measurements. However, both SMC and DMC are likely to explain some aspect of motor control and elevated metabolic power in children with CP.
Despite modest total effects, spasticity itself may not be a large contributor to high metabolic power. In a study by Zaino et al. (2020), children who underwent selective dorsal rhizotomy experienced a 0.9 point reduction in lower body spasticity (approximately a change of one level out of four typically seen) and a 12% reduction in metabolic power, on average. For the same reduction in spasticity, our results indicate that metabolic power should decrease by 2-14% for a child with the median metabolic power. Our estimates capture the values found by Zaino et al. Similar to Zaino's study, our results reflect the total effect of spasticity on metabolic power. Total effects include indirect effects, such as those associated with changes to gait kinematics, so spasticity alone may not be an important determinant of metabolic power. To further support this point, Zaino et al. found changes in metabolic power were not significantly different when compared to a control group, and when we calculate the direct effect of spasticity, it is only a third of its total effect. Therefore, other factors may be responsible for meaningful reductions in metabolic power following spasticity reduction.
Compared to the other variables assessed in this study, strength has smaller causal effects and training strength is not likely to cause meaningful reductions in metabolic power. With a single standard deviation change in strength, our model predicts that metabolic power changes by 1-7% for the average child in our sample. While this may seem like a small reduction, it is not entirely unexpected. As a result of strength training, children with CP often walk faster and for longer distances [24,25]. Since walking faster ultimately increases metabolic power, it would explain why strength has a smaller total effect on metabolic power. Interestingly, strength shares a slight inverted-U-shaped relationship with metabolic power. Weaker individuals may experience increased metabolic power after strength training because they can walk longer and faster (Damiano and Abel, 1998;Eagleton et al., 2004). And stronger individuals may experience decreased metabolic power after strength training because their motor control has improved along with their strength from repetition of movements, enabling them to walk more efficiently [22,23]. Other studies have also speculated that different responses to strength training arise due to changes in walking speed and mechanical efficiency [24,25]. But, although metabolic power may not change drastically, strength training is still important for improving other functional outcomes important to children with CP [54].
We did not normalize or non-dimensionalize metabolic power. We made this choice to simplify our causal model and avoid spurious relationships caused by imperfect normalization. Most non-dimensionalization schemes assume a linear dependence of metabolic power on mass. While this assumption holds true in typically developed populations, it may not apply to children with CP. A study by Plasschaert et al. simulated weight gain during walking in children with and without CP-they noticed that mass normalization was not as effective in removing mass-dependence for the children with CP [58]. We noticed a similar phenomenon in our data, suggesting that mass and metabolic power share a non-linear relationship in children with CP. If this is true, decreases in mass may elicit greater reductions in metabolic power for a child with CP than a typically developed child.
We interpret the effects of walking speed, mass, age, and height to facilitate understanding of our results. Improved walking speed is often an objective of treatments. However, both long term and instantaneous changes in speed come about due to changes in other factors, such as strength, motor coordination or walking mechanics [24,54,59,60]. In our causal model, walking speed was affected by many variables, but only affected GDI [61] and metabolic power [62][63][64]. Since we adjusted for GDI (see Table A in S1 Appendix) and no other variables were affected by walking speed, the total effect of speed is the same as its direct effects. Thus, our results for speed only reflect changes in metabolic power that may occur from walking faster or slower, not from changes in upstream variables. It would be incorrect to compare these results to those found in intervention studies. Both age and height also had relatively large effects on metabolic power but mass was a mediator for both of these relationships. Most of the age and height effects may be explained by changes in mass [65,66]. However, while both mass and speed greatly affect metabolic power, this is the case for all vertebrates and so these variables are unlikely to explain why children with CP have elevated metabolic power compared to typically developing peers.

Limitations
Data reduction enabled us to create a simple causal model, but this also limits our understanding of the specific contributions of the factors that were summarized in the analysis. For example, the GDI is a summary measure of lower body kinematics, so we cannot ascertain if knee flexion or ankle plantarflexion during gait are more important in determining metabolic power. Similarly, SMC and DMC provide summary measures of motor control but quantifying motor control remains challenging. The SMC and DMC measures were available from clinical exams, but only the face validity of SMC and its similarity to other published metrics [67] provide the basis for its use. The power of a DAG is articulating assumed causal relationships in a manner that can be mathematically tested. It is possible that there are direct or indirect causal relationships that we missed in this study, or other variables that may influence responses. In particular, variables related to bony deformities and muscle morphology should be explored further in the future. We did not have data regarding muscle morphology, but it is highly likely to influence muscle strength, kinematics, and metabolic power in children with CP and would likely influence the adjustment sets and effect magnitudes found within this study. It will be important to determine if the addition of such variables would greatly alter the results since the downstream effects of altered muscle morphology may already have been included within the data through differences in strength and kinematics of these children. In this way, the variables included in our causal model may not represent a complete set of variables that cause high metabolic power in CP. However, this DAG did demonstrate good alignment with the data available (i.e., conditional independencies were low) and provides a base for future comparisons or other theories of factors contributing to metabolic power. Future studies should look to improve upon this work.
Causal modeling is a controversial topic. However, controversial does not mean erroneous. Our model, like all models, requires assumptions. An advantage of structural causal models is that the assumptions are made clear and explicit, and the plausibility of those assumptions is tested. So, while our model may not perfectly represent the causes of high metabolic power in CP, by explicitly specifying our causal assumptions, researchers will be able to compare their results and create more comprehensive models in the future. We believe our causal model is a valuable starting point and will lead to experimental research that accepts or rejects hypotheses informed by these results.

Conclusion
Reducing metabolic power can play an important role in reducing or delaying fatigue onset in children with CP and allowing them to be more physically active. Using large-scale data spanning a diverse study sample of widely varying age, GMFCS level, and CP subtype, we sought to understand the causal effects of commonly treated variables (i.e., GDI, SMC, DMC, spasticity, and strength) on metabolic power. Our findings suggest that impaired gait kinematics contribute the most to high metabolic power, followed by SMC and DMC. If the primary goal of treatment is to reduce fatigue, improving motor control and gait mechanics may be the most effective treatment options. Future research should look to perform controlled experiments to test the legitimacy of these findings and possibly investigate other factors that were not reported here.
Supporting information S1 Data. Excel file containing minimal dataset. (XLSX) S1 Appendix. Causal effects of age, height, mass, walking speed, and sex.