Intra- and inter-subject variability of femoral growth plate stresses in typically developing children and children with cerebral palsy

Little is known about the influence of mechanical loading on growth plate stresses and femoral growth. A multi-scale workflow based on musculoskeletal simulations and mechanobiological finite element (FE) analysis can be used to estimate growth plate loading and femoral growth trends. Personalizing the model in this workflow is time-consuming and therefore previous studies included small sample sizes (N < 4) or generic finite element models. The aim of this study was to develop a semi-automated toolbox to perform this workflow and to quantify intra-subject variability in growth plate stresses in 13 typically developing (TD) children and 12 children with cerebral palsy (CP). Additionally, we investigated the influence of the musculoskeletal model and the chosen material properties on the simulation results. Intra-subject variability in growth plate stresses was higher in cerebral palsy than in typically developing children. The highest osteogenic index (OI) was observed in the posterior region in 62% of the TD femurs while in children with CP the lateral region was the most common (50%). A representative reference osteogenic index distribution heatmap generated from data of 26 TD children’s femurs showed a ring shape with low values in the center region and high values at the border of the growth plate. Our simulation results can be used as reference values for further investigations. Furthermore, the code of the developed GP-Tool (“Growth Prediction-Tool”) is freely available on GitHub (https://github.com/WilliKoller/GP-Tool) to enable peers to conduct mechanobiological growth studies with larger sample sizes to improve our understanding of femoral growth and to support clinical decision making in the near future.

1 Introduction

During growth, bones of mammals do not only grow in size but also the geometry adapts according to the predominant loading conditions (Arkin and Katz, 1956; Rauch, 2005; Carter and Beaupré, 2007; Mirtz et al., 2011; Mellon and Tanner, 2012). Altered forces acting on bones can lead to various skeletal deformities. In humans, femoral deformities are common in children with cerebral palsy (CP) but also in children without neurological disorders (Bobroff et al., 1999; Morrell et al., 2002; Beals, 2008). These children are born with typical bony geometry, but in many patients the femoral neck shaft angle (NSA) and anteversion angle (AVA) does not decrease during growth as in typically developing (TD) children (Fabry et al., 1973; Bobroff et al., 1999; Robin et al., 2008). The pathological femoral geometry alters moment arms of muscles (Arnold et al., 2008) which influences muscle recruitment strategy (Kainz et al., 2021b), resulting in altered joint contact forces which are often associated with joint pain (Mackay et al., 2021). Altered movement strategies and joint loading subsequently lead to further progressive disabilities with walking impairments being one of the most profound (Rethlefsen et al., 2017). The altered loading on the growth plate is assumed to be responsible for the pathological femoral growth, yet this assumption remains to be confirmed.

Multi-scale simulations based on musculoskeletal (MSK) simulations and mechanobiological finite element (FE) analysis are used to estimate growth plate loading and femoral growth trends, i.e. the change of the femoral NSA and the AVA due to the forces acting on it (Shefelbine and Carter, 2004; Carriero et al., 2011; Yadav et al., 2016; Yadav et al., 2017; Kainz et al., 2020). First, three-dimensional gait analysis data is collected and MSK simulations are performed to estimate muscle and joint contact forces. Subsequently, these forces are used as loading conditions in a FE model of the femur to quantify the stresses at the growth plate. Based on the estimated shear and hydrostatic stresses at the growth plate the osteogenic index (OI) can be defined, which determines the subject-specific growth rate due to mechanical loading (Carter and Wong, 1988). The mechanical growth rate is added to the biological growth rate which is determined by biological factors, e.g. genetics and nutrients (Stevens et al., 1999). In regions with high OI, bone growth is more likely to occur, whereas in regions with negative OI, growth will be inhibited (Yadav et al., 2016). With these analyses, the development of pathological femoral growth could be identified at an early stage and simple and non-invasive interventions (e.g., orthoses, physical therapy) could potentially normalize bone growth and prevent the necessity of severe corrective osteotomies at a later stage.

Hip joint contact forces (HCF) and subject-specific growth plate geometry are the two main factors influencing growth plate stresses and femoral growth predictions. HCF have been shown to be the main biomarker which determines femoral bone growth (Carriero et al., 2011; Kainz et al., 2020; Yadav et al., 2021). Kainz et al. (2021a) recently showed that the orientation of the HCF in the sagittal plane can differentiate between children with CP who are likely to have typical and pathological femoral growth. HCF are influenced by the subject-specific gait pattern and femoral geometry (Kainz et al., 2020; Kainz et al., 2021b; Modenese et al., 2021). Yadav et al. (2016) showed that the growth plate shape and location affect the stress distribution within the growth plate and therefore change the OI. These studies highlight the importance of accounting for the subject-specific gait pattern and using subject-specific MSK and FE models for mechanobiological growth simulations.

All previous femoral mechanobiological growth studies were based either on a generic FE model of the femur (i.e., generic femoral morphology and growth plate shape and location) (Carriero et al., 2011; Kainz et al., 2020; 2021a) or include very small sample sizes (N < 4) (Yadav et al., 2016; Yadav et al., 2017; Yadav et al., 2021). The reason for this might be the time-consuming generation of personalized FE models. The mechanobiological model requires hexahedral elements aligned with the growth plate. The creation of such a hexahedral mesh from irregular 3D volumes of the femur is challenging and time-consuming. Furthermore, the measured material properties reported in literature (Linde et al., 1985; Rho et al., 1993) for bone vary substantially. Previous studies (Carriero et al., 2011; Yadav et al., 2016; Yadav et al., 2017; Kainz et al., 2020; Kainz et al., 2021a; Yadav et al., 2021) used material property values on the lower bounds of the reported range. Moreover, the generalization of the findings from previous studies was limited due to small sample sizes and the lack of comprehensive reference dataset for TD children. Hence, a reference data set of OI values from TD children obtained with subject-specific MSK and FE models is needed to quantify the typical OI variability within TD children and facilitate an appropriate interpretation of the findings for future femoral growth studies.

This study included the following two clinical aims: 1) quantify and compare the OI in a reference dataset of TD children and children with CP, 2) quantify the intra- and inter-subject variability of the OI. Additionally, our study included the following two technical aims: 3) evaluate the influence of the chosen MSK model, i.e. generic-scaled versus personalized femoral geometry, and 4) evaluate the effect of the chosen material properties on the simulation results. Furthermore, we analyzed and compared the shape, the location and the orientation of the growth plates between groups (CP and TD). Considering that the variability in CP gait is higher compared to TD gait (Kurz et al., 2012; Davies and Kurz, 2013), we hypothesized that the intra-subject variability in the OI is higher in children with CP compared to TD children. In addition, we hypothesized that the chosen MSK model will influence the OI distribution, whereas the material properties will mainly alter the magnitude but not the distribution of the OI. A toolbox based on freely available software, which enables the generation of the personalized FE model and perform the multi-scale simulations in a simple and comprehensive way was developed and used to address our aims. Magnetic resonance imaging (MRI) and three-dimensional gait analysis (3DGA) data of twelve children with CP and thirteen TD children were used to generate subject-specific FE models and quantify the OI based on different input data and model properties.

2 Materials and methods

2.1 Data collection

MRI and 3DGA data including marker trajectories and ground reaction forces of twelve children diagnosed with CP (10.4 ± 3.8 years old, height: 133.6 ± 16.1 cm, mass: 30.1 ± 10.8 kg) and thirteen TD children (10 ± 2.2 years old, height: 144.5 ± 8.5 cm, mass: 36.8 ± 9.5 kg) were analyzed for this study. All participants walked without walking aids and with a self-selected speed. The data of all CP children and three TD children was captured during a previous study (Kainz et al., 2017) while the data of the remaining ten TD children was additionally collected for the purpose of this study. Ethics approval was obtained from the local ethics committees (University of Vienna, reference number 00578). Data collection of the retrospective analyzed data (CP children and three TD children) is described in detail in Kainz et al. (2017).

MRI images of the additionally recorded data (ten TD children) were collected using a 3T magnetic resonance scanner (MAGNETOM Vida, Siemens, Berlin/Munich, Germany) with a T1 vibe sequence with a voxel size of 0.8 × 0.8 × 0.7 mm. 3DGA-data for these ten TD children were captured on the same day as the MRI images using a 12 camera motion capture system (Vicon Motion Systems, Oxford, UK). The used marker set during the motion capturing was based on the Plug-in-Gait marker set (Kadaba et al., 1990; Davis et al., 1991) with additional clusters of three markers on each thigh and shank segment and an additional marker at the 5th metatarsal head of each foot. Simultaneously, ground reaction forces were acquired using five force plates (Kistler Instrumente, Winterthur, Switzerland). All children performed several gait trials with a self-selected walking speed. Marker trajectories were captured, labelled, and filtered (Butterworth 4th order, 6 Hz low-pass filter) in Nexus 2.12.1 (Vicon Motion System, Oxford, United Kingdom).

2.2 Segmentation of MRIs and quantification of geometrical features

3D Slicer 4.13 (Fedorov et al., 2012) was used to segment MRI images. Each femur was split into five parts similar to previous studies (Kainz et al., 2020)—the proximal trabecular bone, the growth plate, the cortical bone of the shaft, the bone marrow and the distal trabecular bone. STL-files of all parts and additionally a file containing the full femur were exported. The STAPLE-Toolbox of Modenese and Renault (2021) was used to identify the femoral head and the epicondyles representing the hip and the knee joint axis using the “GIBOC-Femur” and “GIBOC-Cylinder” algorithms, respectively. If “GIBOC-Cylinder” algorithm failed to fit a cylinder through both epicondyles, “GIBOC-Ellipsoids” algorithm was used to fit ellipsoids through medial and lateral epicondyles. The hip joint center and knee joint axis were required to transform the femur into the OpenSim coordinate system.

The diaphysis of the femur was defined by removing 20% off the top and bottom of the femur. Then, the principal inertia axis of the remaining part was calculated to identify the shaft axis. The neck axis was defined by fitting a least-squares cylinder through surface nodes of the femoral neck. The longitudinal axis of this cylinder was constrained to pass through the femoral head center. The AVA was calculated as the angle between the neck axis and the medial-lateral knee axis obtained from STAPLE-Toolbox (Modenese and Renault, 2021) in the transverse plane. The NSA was computed as the angle between the neck axis and shaft axis in 3D space.

2.3 MSK simulations and loading for FE models

Two models, i.e., generic-scaled and personalized, were created for each participant and used to perform MSK simulations with OpenSim 4.2 to estimate muscle and joint contact forces acting on the femur (Delp et al., 2007; Steele et al., 2012). For the generic-scaled models, the generic ‘gait2392’ OpenSim model (Delp et al., 2007) with locked metatarsophalangeal joints was linearly scaled to fit to the participants’ anthropometry based on the location of surface markers (Kainz et al., 2017). For the personalized model, the Torsion Tool (Veerkamp et al., 2021) was used to modify the femoral geometry in the ‘gait2392’ model to match each child’s NSA and AVA before scaling the model. The maximum isometric muscle forces of all models were scaled based on the ratio of the body mass between the participant’s model and unscaled reference model (Eq. (1)) (van der Krogt et al., 2016; Kainz et al., 2018). In summary, we had two models for each participant which were exactly equivalent except for the femoral geometry and corresponding muscle paths and attachments of muscles.

${\mathit{F}}_{\mathit{s}\mathit{c}\mathit{a}\mathit{l}\mathit{e}\mathit{d}}={\mathit{F}}_{\mathit{g}\mathit{e}\mathit{n}\mathit{e}\mathit{r}\mathit{i}\mathit{c}}*{\left({\mathit{m}}_{\mathit{s}\mathit{c}\mathit{a}\mathit{l}\mathit{e}\mathit{d}}/{\mathit{m}}_{\mathit{g}\mathit{e}\mathit{n}\mathit{e}\mathit{r}\mathit{i}\mathit{c}}\right)}^{\left(\mathbf{2}/\mathbf{3}\right)}(1)$

All models and the corresponding gait analysis data were used to calculate joint angles, joint moments, muscle forces and joint contact forces using inverse kinematics, inverse dynamics, static optimization by minimizing the sum of squared muscle activations and joint reaction load analyses, respectively. Knee and ankle joint markers were only used for scaling and excluded during inverse kinematics. The remaining markers were weighted equally. Maximum marker errors and root-mean-square errors were accepted if less than 4 cm and 2 cm, respectively, as recommended by OpenSim’s best practice recommendations (Hicks et al., 2015). Additional analyses were performed to identify muscle attachments on the femur and obtain the effective directions of muscle forces (van Arkel et al., 2013). The mean waveform of the resultant HCF from all trials was calculated and the trial with the lowest root mean square difference to the mean waveform was selected as a representative loading condition. Similar to previous studies (Yadav et al., 2016; Kainz et al., 2020) nine load instances were selected based on the HCF peaks and the valley in-between during the stance phase. The HCF and muscle forces acting on the femur during the nine load instances were used as loading conditions for FE analysis.

2.4 Creation of hexahedral mesh

The mechanobiological model requires hexahedral elements aligned within the growth plate stacked in several layers to define transition zones and enable progressive growth simulations for each layer of the growth plate. We developed the GP-Tool to automatically create hexahedral meshes based on the subject-specific femoral geometry. The STL-files obtained from the segmented femurs were used as input for the GP-Tool. A visual overview of the steps to create a mesh with the GP-Tool is shown in Figure 1.

FIGURE 1

FIGURE 1. Visual description of the steps performed by the semi-automated toolbox with (A) the segmentation of the MRI data as input. Subfigure (B) shows the 3D geometry of the femur in the MRI reference system. In subfigure (C) the geometry was rotated to have the growth plate parallel to the XY plane, the growth plate was removed, the resulting gap closed by moving the proximal part of the femoral head and the edges smoothed. Sculpt tool of Coreform Cubit was used to create a hexahedral mesh of this geometry. Subfigure (D) shows the final mesh with ten laysers of elements aligned with the growth plate.

Using MATLAB’s (MathWorks, Natick, MA, United States) principal component analysis (“pca”-function) the main orientation of the growth plate was determined. All parts of the femur were rotated so that growth plates’ main orientation was parallel to the XY plane. The growth plate itself was removed and the part above the growth plate of the proximal trabecular bone was positioned on the part which is below the growth plate. Smoothing of the intersection region was performed with MeshLab (Cignoni et al., 2008). These steps were necessary to create a continuous mesh to add perfectly aligned hexahedral elements in the growth plate later. Sculpt tool of Coreform Cubit (Coreform, Utah, United States) was used to create a hexahedral mesh with an element size of approximately 1.5 mm. A mesh convergence study was conducted based on three femurs to ensure that the results are not influenced by the number of elements in the mesh (see Supplementary Material). The part above the growth plate was moved to its original position and ten layers with equal height were added and presented the growth plate. Finally, the mesh was optimized to no longer include elements with negative Jacobians. This procedure was performed for each femur of all participants (N = 50) and resulted in meshes of approximately 150.000 nodes and 140.000 elements varying due to different femur sizes.

2.5 Finite element analysis

A FE model was created for each femur based on the subject-specific hexahedral mesh and the loading conditions obtained from the MSK simulations. All models were fixed at the femoral epicondyles and HCF and muscle forces were applied as nodal forces for the nine load instances.

For each load instance the HCF was distributed to the closest 100 surface nodes (approximately 2.25 cm²) in the direction of the HCF orientation. For each muscle a node was identified which was the closest node to the muscle attachment obtained from the OpenSim simulations (van Arkel et al., 2013). Due to discrepancies in the geometry, e.g., bending of the shaft, between OpenSim’s femur and the participants’ femur derived from MRI an algorithm was used to ensure that the defined muscle attachment was on the same side (anterior/posterior or medial/lateral) of the femur. Nodal forces were applied to this node in x/y/z directions according to the muscle directions obtained from the additional muscle analysis (van Arkel et al., 2013) in order that the resulting force was equal to the muscle force estimated by the static optimization algorithm.

The FE model was duplicated and two different sets of linear elastic materials with Young’s modulus and Poisson ratio described in Table 1 were assigned to the different parts of the femur. The chosen values for material properties were based on literature and previously used values in mechanobiological growth studies (Linde et al., 1985; Rho et al., 1993; Carriero et al., 2011; Yadav et al., 2016; Kainz et al., 2020). A transition zone of three layers (out of the ten layers within the growth plate) between trabecular bone and the growth plate was modeled with linearly decreasing Young’s modulus from the trabecular bone to the growth plate to represent the mineralizing bone tissue (Kainz et al., 2020). FEBio 3 (Maas et al., 2012) was used for FE simulations and to calculate principal stresses.

TABLE 1

TABLE 1. Sets of material properties used for the different sections of the femur.

2.6 Osteogenic index calculation

The growth rate due to mechanical loading was estimated as the OI (Stevens et al., 1999). Considering that bone mineralizes on the diaphyseal side, the OI was determined for each element of the most distal layer within the growth plate (excluding transition zone). Similar to previous studies (Shefelbine and Carter, 2004; Carriero et al., 2011; Yadav et al., 2016; Yadav et al., 2017; Kainz et al., 2020; Kainz et al., 2021a; Yadav et al., 2021) the OI was estimated using the obtained principal stresses (σ₁, σ₂, σ₃) from the FE analysis (Eq. (2)). Hydrostatic stresses (${\sigma }_{H_i}$) and octahedral shear stresses (${\sigma }_{S_i}$) were calculated for each load instance (i = 1…9) (Equations (3) and (4) and constants a and b were chosen with values of 0.02 MPa⁻¹month⁻¹ and 0.01 MPa⁻¹month⁻¹, respectively, to have a ratio b/a = 0.5 similar to previous studies) (Shefelbine and Carter, 2004; Carriero et al., 2011; Yadav et al., 2016; Yadav et al., 2017; Yadav et al., 2021). The unit of the OI is therefore month⁻¹. Outliers identified as OI values that were more than three scaled median absolute deviations away from the median (MATLAB’s “isoutlier” function) were replaced by the corresponding minimum or maximum values. These identified outliers were located at the outer edges of the growth plate, which might be due to the lower quality of the hexahedral elements in this region compared to the rest of the growth plate. Therefore, this step was necessary to ensure smoother OI distribution for subsequent comparisons.

$\mathit{O}\mathit{I}=\mathit{a}\ast \mathit{max}\left({\mathit{\sigma }}_{{\mathit{S}}_{\mathit{i}}}\right)+\mathit{b}\ast \mathit{min}\left({\mathit{\sigma }}_{{\mathit{H}}_{\mathit{i}}}\right)(2)$

${\mathit{\sigma }}_{{\mathit{S}}_{\mathit{i}}}=\frac{\sqrt{{\left({\mathit{\sigma }}_{\mathbf{1}}-{\mathit{\sigma }}_{\mathbf{2}}\right)}^{\mathbf{2}}+{\left({\mathit{\sigma }}_{\mathbf{2}}-{\mathit{\sigma }}_{\mathbf{3}}\right)}^{\mathbf{2}}+{\left({\mathit{\sigma }}_{\mathbf{3}}-{\mathit{\sigma }}_{\mathbf{1}}\right)}^{\mathbf{2}}}}{\mathbf{3}}(3)$

${\mathit{\sigma }}_{{\mathit{H}}_{\mathit{i}}}=\frac{{\mathit{\sigma }}_{\mathbf{1}}+{\mathit{\sigma }}_{\mathbf{2}}+{\mathit{\sigma }}_{\mathbf{3}}}{\mathbf{3}}(4)$

2.7 Quantification of the osteogenic index

To enable a comparison of OI between growth plates with different shapes, the OI values were projected on the transverse plane according to the elements’ locations and interpolated to a squared grid. A blue to red color scheme was used to visualize OI values representing minimum to maximum OI values, respectively. This resulted in heatmaps of equal size for all growth plates revealing the OI distribution (anterior/posterior and medial/lateral) based on the mechanical loading condition. Subsequently, the squared grid was divided into five regions representing the center, anterior, posterior, lateral and medial sections of the growth plate. The center was defined as a circle with a diameter of 50% of the growth plate’s width and the other regions according to their corresponding anatomical sides of the growth plate separated by the square’s diagonals (Figure 2). Consecutively, the region with the highest and lowest mean value was identified which indicates the area with maximally promoted and inhibited growth, respectively. Since these heatmaps were scaled to each analysis’ minimum and maximum values the magnitude of the OI was neglected. Therefore, we additionally evaluated the range of the magnitude and the mean and median value of the OI.

FIGURE 2

FIGURE 2. An example of a projection of OI values on the transverse plane showing the division into five regions. The OI distribution is represented using a blue to red color scheme representing low and high values, respectively.

2.8 Osteogenic index reference dataset

A representative reference OI distribution heatmap was created by calculating mean OI values across participants for the TD and the CP group. One can imagine this as stacking up all heatmaps. Due to the different shapes of the growth plates, on the outer edges only few OI values were available. To avoid incorrect values at these locations, mean values were only calculated for locations if data of more than half of the femurs, i.e. more than 12 and 13 for CP and TD, respectively, were available.

2.9 Intra-subject variability

Two main metrics were used to quantify the intra-subject variability of the OI for each participant. Firstly, the regions with the highest and lowest mean value were compared between each participant’s left and right femur. Secondly, to quantify differences in shape and distribution of the OI (e.g. ring shape vs. linear-gradient), an image comparison between the heatmaps of the left and right femur was performed using OpenCV’s template matching (Bradski, 2000). The “TM_CCOEFF_NORMED”-method was used which computes one resultant value between 0 and 1 that represents the probability that images are equal. For a more comprehensible metric, the result was subtracted from 1 so that high values indicate high variability. Additionally, to account for the magnitude of the OI, the difference of the mean values and difference of the ranges between the left and right OIs were calculated. Subsequently, independent t-tests were performed to identify differences of the intra-subject variability between the CP and TD group.

2.10 Inter-subject variability

To assess the inter-subject variability between groups, we analyzed the mean occurring OI values in each region. Firstly, the mean value of each region was normalized to the range of the individual OI. Then, these normalized values for each region were compared between CP and TD group with independent t-tests to evaluate differences in general OI distribution. To identify statistically significant differences in the variability, the standard deviations of each region were compared between groups with a Levene’s test.

Additionally, we compared each OI heatmap with all OI heatmaps of the same group, i.e. each OI of a TD child was compared to all other OIs of TD children, with the before mentioned image comparison method. Comparisons of equal heatmaps (e.g. “TD01 left side” vs. “TD01 left side”) were neglected. Subsequently, the mean variability to other OIs of the group was calculated. This measure indicates the difference of shape and distribution of a single OI within its group. Next, the inter-subject variability within each group was compared between TD children and children with CP. Again, to account for the magnitude of the OI, the mean values and the ranges of the OIs were compared between groups with a two-sampled t-test.

2.11 Influence of MSK model and FE material properties

The sum of all muscle forces acting on the femur as well as the magnitude of the HCF at the first and the second peak were compared between the MSK simulations with generic-scaled and personalized models. Within each group (TD and CP), all parameters were tested for significant differences with two-sampled t-tests. Between groups, muscle forces and HCF of the personalized MSK model were compared with independent t-tests.

To answer our technical research questions on the influence of the MSK model and the material properties on the OI, we performed the same comparisons and statistical tests as described above for the simulation results obtained with altered loading conditions and altered material properties, separately. Subsequently, we identified differences in outcomes due to different loading conditions and material properties. Additionally, we evaluated changes of the regions with the highest mean value compared to the base analysis with “hard” material properties and personalized MSK models.

2.12 Geometrical analysis

To identify possible reasons for observed differences of the OI between the CP and the TD group, the shape, the location and the orientation of the growth plates were analyzed. The shape was quantified by fitting a sphere through nodes in the middle of the growth plate. If the sphere’s midpoint was distal of the growth plate the shape was convex, otherwise it was concave. The radii of these spheres were compared between groups. The location of the growth plate was determined by calculating the distance between the intersection point of the neck axis and the femur’s proximal surface and the mean location of all nodes within the growth plate. The neck length was quantified as the femoral head radius plus the distance between the femoral head center and the point where neck and shaft axis have the minimal distance to each other. The location of the growth plate was normalized to the length of the femoral neck. To quantify the orientation of the growth plate the angles between the normal vector to the growth plate (obtained during mesh creation) and 1) the femur coordinate system and 2) the orientation of the mean and maximum HCF during stance phase in the anatomical planes were calculated. Additionally, angles between the neck axis and 1) the orientation of the mean and maximum HCF in the anatomical planes and 2) the normal vector to the growth plate were computed.

3 Results

3.1 Participant exclusion and FE analysis metrics

In total, 6 out of 50 (12%) FE analysis with “hard” material properties and the personalized MSK model did not converge and therefore four participants (all children with CP; in two participants, left and right side did not converge) had to be excluded from our study because intra-subject comparisons were impossible. Data of these participants was not used in any analysis. Therefore, all presented results include data of 8 children with CP and 13 TD children, i.e. 42 femurs in total. In the analysis sets with the generic-scaled MSK model and “soft” material properties 2 (4.75%) and 4 (9.5%) out of the remaining 42 FE analysis did not converge, respectively.

3.2 Osteogenic index reference dataset

The representative reference OI distribution heatmap generated from data of 26 TD children’s femurs showed a ring shape with low values in the center region and high values at the border of the growth plate (Figure 3A). The descending order of the regions according to their highest mean was posterior, medial, lateral, anterior and center. The reference OI for TD children was distributed between −0.024 months⁻¹ and 0.013 months⁻¹ with a mean value of −0.003 months⁻¹.

FIGURE 3

FIGURE 3. Average growth plate shape and OI distribution in TD children (A) and children with CP (B). The OI distribution is represented using a blue to red color scheme representing low to high values, respectively.

The reference OI distribution generated from data of 16 femurs of children with CP had a different shape with highest values on the lateral side of the growth plate (Figure 3B). An indication of a ring shape distribution, similar to the TD children, was observed. The regions with highest mean values in descending order were lateral, posterior, anterior and medial followed by the center. The OI values ranged from −0.011 months⁻¹ to 0.022 months⁻¹ with a mean of 0.003 months⁻¹.

3.3 Intra-subject variability

When comparing the left with the right growth plate of each participant, the highest mean OI value occurred at different regions in 30.8% (n = 4) of TD children and in 62.5% (n = 5) of children with CP (Figure 4). The lowest mean OI value was observed at the center region in all TD children and most children with CP. In one child with CP the lowest mean OI value occurred in the medial region and in another child with CP in the posterior region.

FIGURE 4

FIGURE 4. Intra-subject comparison of region with highest mean OI value between left and right growth plate for TD children (left) and children with CP (right).

The differences of the mean OI value as well as the range of OI values between the left and right side were not significantly different between both groups. The image comparison of the heatmaps to quantify the variability of the shape and OI distribution between left and right growth plate showed a significantly higher (p < 0.01) intra-subject variability in children with CP compared to TD children (Figure 5).

FIGURE 5

FIGURE 5. Comparison of the intra-subject variability of the magnitude of OI, i.e. mean value and range, and the variability assessed with template matching of the heatmaps of the left and right OI between children with CP and TD children. Higher variability values indicate a higher difference in shape and distribution, a value of 0 indicates equal OI shape and distribution. Significant differences were quantified with independent t-tests.

3.4 Inter-subject variability

In 61.5% (n = 16) of the femurs of the TD group the highest mean value of the OI was observed in the posterior region of the growth plate followed by the medial region (27%; n = 7). In children with CP the lateral (50%; n = 8) followed by the posterior (31%; n = 5) region were the areas with the most occurred highest mean OI values (Figure 6 left).

FIGURE 6

FIGURE 6. The left diagram shows how often the highest mean value was observed in a specific region. In the right diagram the colored bars represent the mean value of each region normalized to the range of the individual OI. Significant differences were quantified with independent t-tests. Vertical error bars indicate the standard deviation between participants within groups.

The mean values of the posterior and medial regions were significantly higher in TD children compared to children with CP whereas significantly lower values were observed in the center region for the TD children. Levene’s test showed significant differences of the variances (p < 0.05) in the posterior and medial region between groups indicating that in children with CP the variability is higher compared to TD children (Figure 6 right).

The range of the OI as well as the mean OI differed significantly (p < 0.01) between the TD and CP group. In TD children, the range of the OI was higher and lower mean values were observed. Comparing each OI heatmap with all other OI heatmaps of its group revealed significantly higher (p < 0.001) inter-subject variability in the CP compared to the TD group (Figure 7).

FIGURE 7

FIGURE 7. Comparison of the inter-subject variability of the magnitude of OI, i.e. mean value and range, and the variability assessed with template matching of the heatmaps. Each OI was matched with all other OIs of its group and the mean variability compared between groups. Higher variability values indicate a higher difference in shape and distribution, a value of 0 indicates equality. Significant differences were quantified with independent t-tests.

3.5 Geometrical analysis

The femurs of TD children were significantly (p < 0.001) larger compared to those of children with CP. No significant differences were observed for the neck length and the growth plate location between groups normalized to the femur’s length and its neck, respectively. A wide range of AVA and NSA was present in the analyzed dataset but without any statistical differences between the CP and TD group. The shape of the growth plate did not differ between groups, all had a convex shape with radii between 10mm and 35 mm (Figure 8).

FIGURE 8

FIGURE 8. Comparison of (A) the femur size, (B) the neck length, (C) the growth plate location, (D) the shape of the growth plate, (E) the AVA and (F) the NSA between children with CP and TD children. Growth plate location was assessed as the distance between the intersection point of the neck axis and the proximal surface and the middle of the growth plate. The shape of the growth plate was assessed by fitting a sphere into the nodes of the middle layer of the growth plate. Significant differences were quantified with independent t-tests.

The orientation of the growth plate was similar in relation to the femur’s coordinate system between children with CP and TD children. The angles between the orientation of the maximum and mean occurring HCF during stance phase and the vector normal to the growth plate did not differ significantly between groups (Figure 9).

FIGURE 9

FIGURE 9. Comparison of the growth plate orientation (normal vector) between the CP and TD group in relation to the femur’s coordinate system (top), the maximum HCF (middle) and the mean HCF during stance phase (bottom) in transverse, frontal and sagittal plane. Significant differences were quantified with independent t-tests.

No significant difference was observed between the CP and TD group for the angle between the neck axis and the normal vector to the growth plate. The angle between the femoral neck axis and the orientation of the maximum (p < 0.05) and mean (p < 0.001) HCF was significantly different between CP and TD children in transverse plane. A statistically significant difference between CP and TD children was also observed in the frontal and sagittal plane for the angle between the neck axis and mean HCF during stance phase (Figure 10).

FIGURE 10

FIGURE 10. Comparison of the neck axis orientation between the CP and TD group in relation to the growth plates normal vector (top), the orientation of the maximum HCF (middle) and the mean HCF during stance phase (bottom) in transverse, frontal and sagittal plane. Significant differences were quantified with independent t-tests.

3.6 Influence of MSK model

The use of a generic-scaled instead of personalized MSK model changed the loading condition, i.e. HCF and muscle forces on the participants’ femurs (Figure 11). At the first peak of the HCF, the magnitude of the HCF and the sum of all muscles acting on the femur were significantly lower in both groups when using a generic-scaled model. In TD participants, the magnitude of the HCF and the muscle forces at the second peak of HCF were significantly lower with the generic-scaled compared to the personalized model. Independent t-tests showed significantly higher values of HCF peak as well as the sum of all muscle forces acting on the femur at both HCF peaks for TD compared to CP children. Two out of the 42 FE analysis did not converge when using the generic-scaled loading conditions. Against the general observation that HCF and muscle forces were lower when using a generic-scaled MSK model, in these two analyses the maximum HCF magnitude increased which led to non-convergence in the FE analysis.

FIGURE 11

FIGURE 11. Comparison of the HCF and the muscle forces acting on the participants’ femurs between the generic-scaled and the personalized MSK model within the CP and TD group as well as between groups for the personalized MSK model. Significant differences were quantified with dependent and independent t-tests within groups and between groups, respectively.

The region with the highest mean OI value changed due to the MSK model that was used in 2 (14.3%) and 5 (19.2%) participants of the CP and TD group, respectively (Figure 12). The intra-subject comparison between left and right OI showed that in 30.8% (n = 4) of TD children and in 71.4% (n = 5 out of 7) of children with CP the highest mean OI value occurred in different regions of the growth plate. No statistically significant changes were observed for the differences of the mean OI, the range of the OI and the OI variability between the left and right side due to the MSK model that was used. The range of OI did not significantly differ between the CP and TD group when using the loads from the generic-scaled MSK models, whereas it was significantly different with personalized MSK models. Additional plots, similar to Figure 5; Figure 6 but based on the generic-scaled MSK model can be found in the Supplementary Material.

FIGURE 12

FIGURE 12. Comparison of region with highest mean value between the generic-scaled and personalized MSK model for TD children (left) and children with CP (right).

3.7 Influence of FE material properties

For four children with CP the FE analysis did not converge with the “soft” material properties. This reduced the number of intra-subject comparisons for CP children to 5 participants. The region with the highest mean OI value changed due to the material properties in 2 (16.7%) and 5 (19.2%) participants of the CP and TD group, respectively (Figure 13). The intra-subject comparison between left and right OI showed that in 38.5% (n = 4) of TD children and in 60% (n = 3 out of 5) of children with CP the highest mean OI value occurs in different regions of the growth plate. With “soft” material properties, there was no significant difference of the intra-subject variability of the shape and distribution between the left and right side. Inter-subject comparison between the CP and TD group with “soft” material properties showed no significant differences in mean values for all regions. The mean OI was not significantly different between groups whereas this was the case for “hard” material properties. Additional plots, similar to Figure 5; Figure 6 but based on the “soft” material can be found in the Supplementary Material.

FIGURE 13

FIGURE 13. Comparison of region with highest mean value between FE analysis with “hard” and “soft” material properties for TD children (left) and children with CP (right).

4 Discussion

The aim of this study was to quantify and compare the OI in the proximal growth plate of the femur in a comparably large cohort of TD children and children with CP with subject-specific loading, femoral morphology as well as growth plate shape and orientation. The GP-Tool based on freely available software was developed, which enables one to create subject-specific FE models and perform the multi-scale simulations in a simple and comprehensive way. The GP-Tool has been shared with the community (https://github.com/WilliKoller/GP-Tool) to enable peers to conduct mechanobiological growth studies with larger sample sizes and enhance our understanding of femoral bone growth. In this study, the analyzed dataset included 26 femurs of TD children and 16 femurs of children with CP and empowered us to create a reference OI as well as to evaluate the intra- and inter-subject variability. The within-group variability in range of OI, mean OI, OI shape and distribution were significantly different between groups, suggesting different bone growth patterns may be present in CP and TD children.

The reference OI for TD children showed a clear ring shape with high values on the outside and low values in the center which indicates more growth on the outside compared to the inner part of the growth plate. This could support the finding that with time and under normal loading conditions, the growth plate changes its shape from a spherical to a flat disc (Kandzierski et al., 2012).

In children with CP the reference OI clearly differed from those of TD children with the maximum values occurring in the lateral region. Considering the high variability of femoral morphology and femoral loading patterns in children with CP, the presented reference OI does not reflect the whole population of children with CP and therefore should be interpreted with caution. Some children with CP had OIs that were similar to the reference OI of TD children while others were very different.

The inter-subject analyses showed that in most femurs of TD children, the region with the highest mean OI values were the posterior followed by the medial region. In both regions, we found significantly higher OI values within the TD compared to the CP group. Most previous studies based on loading conditions estimated from gait of TD children also found a ring-shaped OI distribution with the highest values on the posterior side of the growth plate (Yadav et al., 2017; Kainz et al., 2020). Carriero et al. (2011) reported a linear-gradient OI distribution with the highest values on the lateral side for a healthy child, which was contrary to the OI of our TD children but similar to the OI of our participants with CP. Carriero et al. (2011) however, conducted their simulations based on a model with a simplified growth plate, i.e. flat disc. The low participant number and simplified growth plate (Yadav et al., 2016) might explain the contrary results from Carriero et al. (2011) compared to the more recent investigations.

The analysis of the intra-subject variability showed that in most TD children the region with the highest mean value was equal between the left and right side while in the CP group this was only the case for approximately a third of the participants. The variability of the shape and distribution of the OI also revealed a significantly higher intra-subject variability in children with CP compared to TD children which could lead to asymmetric femoral growth.

In the analyzed dataset, the femurs of children with CP were significantly smaller compared to the TD group but the morphology in terms of neck length, AVA, NSA, growth plate location and growth plate shape was similar between groups. All morphological parameters were within the range of TD values reported in a study based on 508 participants (Szuper et al., 2015).

The use of generic-scaled instead of image-based MSK models significantly changed the femoral loading conditions, i.e. muscle forces and HCF, but the general findings remained similar. The only statistical differences between analysis based on the generic-scaled and personalized MSK model were observed for the range of the OI, whereas the generic-scaled MSK model led to fewer differences between the CP and TD group.

Four analyses did not converge when using material properties with the compliant material, i.e., low Young’s modulus. In linear elastic materials, which were used in this study, lower stiffness leads to larger displacements. When elements are deformed resulting into inverted elements, the FE solver cannot continue the analysis. In theory, the stress within a material only depends on the force and the cross-section of the geometry and therefore, stresses within the growth plate should not be influenced by the used material properties. However, in FE analysis the material properties are used to initially calculate the global stiffness matrix. Then, displacements are estimated and subsequently used to calculate the stresses (Bathe, 2006, 149ff). Hence, changes of material properties have an impact on observed stresses within the growth plate. However, observed changes in our analyses were minimal with the only statistical difference in the mean OI value reducing the difference between CP and TD group.

Four participants had to be excluded because in total, 6 out of 50 FE analysis did not converge with stiff material properties. The reasons for this were mainly tremendously high HCFs (∼4000 N, above 10 times body weight) or very thin cortical bone. Thin cortical bone, i.e., only two or three element layers on the femoral shaft (bone marrow has low Young’s modulus and can be neglected), is unlikely to support the full MSK loading. In future studies, attention has to be paid to the segmentation of the cortical bone. Unrealistic and insufficient cortical thickness and/or unrealistic high HCF might not converge and therefore make the estimation of the OI impossible.

We only estimated growth plate stresses and quantified the OI but did not simulate femoral growth. Femoral growth can be predicted with an additional FE analysis including a method to estimate the growth direction. To do this, different methods have been suggested to define the growth direction, i.e. femoral neck deflection direction or principal stress direction (Yadav et al., 2016). None of these methods have been validated based on longitudinal data and therefore we decided not to include the additional analysis in our study.

In children with CP muscle morphology and motor control are often altered compared to TD children, i.e. shorter muscle belly, muscle spasticity, decreased muscle volume and strength (Hanssen et al., 2021; Peeters et al., 2023). Assessment of individuals’ muscle properties, e.g. maximum isometric muscle forces, with clinical tests is challenging in children with CP (Kainz et al., 2018). We accounted for the individuals’ bony geometry which has been shown to have a larger impact on hip joint contact forces than the personalized motor control (Kainz et al., 2021b) Nevertheless, future studies could include ultrasound measurements to personalize muscle parameters (Schless et al., 2018), include a spasticity model (Falisse et al., 2018; Falisse et al., 2020) and account for the patient-specific motor control (Veerkamp et al., 2019) to further improve the accuracy of the MSK simulations.

Our study included the following limitations. First, bone is a highly complex structure with anisotropic material properties. In this study, the material of the bone was simplified to have linear elastic properties. Previous studies modelled the bone similarly (Shefelbine and Carter, 2004; Carriero et al., 2011; Yadav et al., 2016; Kainz et al., 2020; Yadav et al., 2021). FEBio supports complex materials, therefore, the developed GP-Tool can be modified to include more realistic materials for future studies. Second, during the automated meshing, a small number of elements with low quality metrics (at the very outer border) of the growth plate were removed. This was necessary to avoid inverted elements under load to ensure successful FE analysis. Third, the MSK models included only the personalized femoral morphology whereas the geometry of the tibia was neglected. Changes of bony geometry alters moment arms of muscles and subsequently has an impact on the estimated muscle forces and HCF. However, tibial torsion only affects a limited number of muscles and therefore we assume the impact on HCF is negligible. Studies which investigate the impact of tibial torsion on HCF should be carried out to verify our assumption. Fourth, we used a static optimization algorithm to compute muscle activities. Approaches including electromyography (EMG) data might better represent the individuals’ motor control (Pizzolato et al., 2015). Due to lack of EMG data in all our participants we could not perform EMG-informed analysis and therefore decided to use static optimization. Femoral geometry has a bigger impact on HCF than including EMG data in the optimization (Kainz et al., 2021b) and therefore we would expect to get similar OI with an EMG-informed approach.

The findings of this paper can help to identify possible pathological loading on the growth plate which might result in atypical femoral bone growth. If pathological loading is noticed at an early stage, subject-specific therapies like real-time biofeedback training could be used to alter joint loads (Comellas and Shefelbine, 2022; Uhlrich et al., 2022) and therefore change growth plate stresses. However, further studies have to be carried out to identify how different walking pattern and motor control strategies influence growth plate stresses.

To conclude, we developed the GP-Tool to create subject-specific FE models and estimate the OI in a semi-automatic way. The code of the GP-Tool has been made freely available to the research community as well as an executable file for users without a MATLAB license. On a standard notebook, the complete process of creating a subject-specific mesh from STL files, setting up and running the FE analysis takes around 30 min. We hope the GP-Tool encourages peers to conduct similar studies with large sample sizes to improve our general insights in femoral bone growth and to support clinical decision making in the near future.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving human participants were reviewed and approved by Ethics commitee of the University of Vienna and Children’s Health Queensland Human Research Ethics Committee. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

Author contributions

WK, AB, and HK contributed to conception and design of the study. WK and HK collected and processed data. WK developed the GP-Tool which was used to perform OI analysis. HK, AB, and BG supported the development with ideas about implementation, testing and feedback. WK wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This research was partly funded by a PhD grant from the PhaNuSpo Vienna Doctoral School.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2023.1140527/full#supplementary-material

References

Arkin, A. M., and Katz, J. F. (1956). The effects of pressure on epiphyseal growth: The mechanism of plasticity of growing bone. JBJS 38, 1056–1076. doi:10.2106/00004623-195638050-00009

PubMed Abstract | CrossRef Full Text | Google Scholar

Arnold, A. S., Komallu, A. V., and Delp, S. L. (2008). Internal rotation gait: A compensatory mechanism to restore abduction capacity decreased by bone deformity? Dev. Med. Child Neurol. 39, 40–44. doi:10.1111/j.1469-8749.1997.tb08202.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bathe, K. J. (2006). Finite element procedures. Prentice-Hall. Available at: https://books.google.at/books?id=rWvefGICfO8C.

Google Scholar

Beals, R. K. (2008). Developmental changes in the femur and acetabulum in spastic paraplegia and diplegia. Dev. Med. Child Neurol. 11, 303–313. doi:10.1111/j.1469-8749.1969.tb01437.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bobroff, E. D., Chambers, H. G., Sartoris, D. J., Wyatt, M. P., and Sutherland, D. H. (1999). Femoral anteversion and neck-shaft angle in children with cerebral palsy. Clin. Orthop. Relat. Res. 364, 194–204. doi:10.1097/00003086-199907000-00025

PubMed Abstract | CrossRef Full Text | Google Scholar

Bradski, G. (2000). The OpenCV library. Dr. Dobb’s Journal of Software Tools.

Google Scholar

Carriero, A., Jonkers, I., and Shefelbine, S. J. (2011). Mechanobiological prediction of proximal femoral deformities in children with cerebral palsy. Comput. Methods Biomech. Biomed. Eng. 14, 253–262. doi:10.1080/10255841003682505

PubMed Abstract | CrossRef Full Text | Google Scholar

Carter, D. R., and Beaupré, G. S. (2007). Skeletal function and form: Mechanobiology of skeletal development, aging, and regeneration. Cambridge University Press.

Google Scholar

Carter, D. R., and Wong, M. (1988). The role of mechanical loading histories in the development of diarthrodial joints. J. Orthop. Res. 6, 804–816. doi:10.1002/jor.1100060604

PubMed Abstract | CrossRef Full Text | Google Scholar

Cignoni, P., Callieri, M., Corsini, M., Dellepiane, M., Ganovelli, F., and Ranzuglia, G. (2008). “MeshLab: An open-source mesh processing tool,” in Eurographics italian chapter conference. Editors V. Scarano, R. D. Chiara, and U. Erra (The Eurographics Association). doi:10.2312/LocalChapterEvents/ItalChap/ItalianChapConf2008/129-136

CrossRef Full Text | Google Scholar

Comellas, E., and Shefelbine, S. J. (2022). The role of computational models in mechanobiology of growing bone. Front. Bioeng. Biotechnol. 10, 973788. doi:10.3389/fbioe.2022.973788

PubMed Abstract | CrossRef Full Text | Google Scholar

Davies, B. L., and Kurz, M. J. (2013). Children with cerebral palsy have greater stochastic features present in the variability of their gait kinematics. Res. Dev. Disabil. 34, 3648–3653. doi:10.1016/j.ridd.2013.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Davis, R. B., Õunpuu, S., Tyburski, D., and Gage, J. R. (1991). A gait analysis data collection and reduction technique. Hum. Mov. Sci. 10, 575–587. doi:10.1016/0167-9457(91)90046-Z

CrossRef Full Text | Google Scholar

Delp, S. L., Anderson, F. C., Arnold, A. S., Loan, P., Habib, A., John, C. T., et al. (2007). OpenSim: Open-Source software to create and analyze dynamic simulations of movement. IEEE Trans. Biomed. Eng. 54, 1940–1950. doi:10.1109/TBME.2007.901024

PubMed Abstract | CrossRef Full Text | Google Scholar

Fabry, G., Macewen, G. D., and Shands, A. R. (1973). Torsion of the femur: A follow-up study in normal and abnormal conditions. J. Bone Jt. Surg. 55, 1726–1738. doi:10.2106/00004623-197355080-00017

PubMed Abstract | CrossRef Full Text | Google Scholar

Falisse, A., Bar-On, L., Desloovere, K., Jonkers, I., and De Groote, F. (2018). A spasticity model based on feedback from muscle force explains muscle activity during passive stretches and gait in children with cerebral palsy. PLoS ONE 13, e0208811. doi:10.1371/journal.pone.0208811

PubMed Abstract | CrossRef Full Text | Google Scholar

Falisse, A., Pitto, L., Kainz, H., Hoang, H., Wesseling, M., Van Rossom, S., et al. (2020). Physics-based simulations to predict the differential effects of motor control and musculoskeletal deficits on gait dysfunction in cerebral palsy: A retrospective case study. Front. Hum. Neurosci. 14, 40. doi:10.3389/fnhum.2020.00040

PubMed Abstract | CrossRef Full Text | Google Scholar

Fedorov, A., Beichel, R., Kalpathy-Cramer, J., Finet, J., Fillion-Robin, J.-C., Pujol, S., et al. (2012). 3D slicer as an image computing platform for the quantitative imaging network. Magn. Reson Imaging 30, 1323–1341. doi:10.1016/j.mri.2012.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanssen, B., Peeters, N., Vandekerckhove, I., De Beukelaer, N., Bar-On, L., Molenaers, G., et al. (2021). The contribution of decreased muscle size to muscle weakness in children with spastic cerebral palsy. Front. Neurol. 12, 692582. doi:10.3389/fneur.2021.692582

PubMed Abstract | CrossRef Full Text | Google Scholar

Hicks, J. L., Uchida, T. K., Seth, A., Rajagopal, A., and Delp, S. L. (2015). Is my model good enough? Best practices for verification and validation of musculoskeletal models and simulations of movement. J. Biomech. Eng. 137, 020905. doi:10.1115/1.4029304

PubMed Abstract | CrossRef Full Text | Google Scholar

Kadaba, M. P., Ramakrishnan, H. K., and Wootten, M. E. (1990). Measurement of lower extremity kinematics during level walking. J. Orthop. Res. 8, 383–392. doi:10.1002/jor.1100080310

PubMed Abstract | CrossRef Full Text | Google Scholar

Kainz, H., Hoang, H. X., Stockton, C., Boyd, R. R., Lloyd, D. G., and Carty, C. P. (2017). Accuracy and reliability of marker-based approaches to scale the pelvis, thigh, and shank segments in musculoskeletal models. J. Appl. Biomech. 33, 354–360. doi:10.1123/jab.2016-0282

PubMed Abstract | CrossRef Full Text | Google Scholar

Kainz, H., Goudriaan, M., Falisse, A., Huenaerts, C., Desloovere, K., De Groote, F., et al. (2018). The influence of maximum isometric muscle force scaling on estimated muscle forces from musculoskeletal models of children with cerebral palsy. Gait Posture 65, 213–220. doi:10.1016/j.gaitpost.2018.07.172

PubMed Abstract | CrossRef Full Text | Google Scholar

Kainz, H., Killen, B. A., Wesseling, M., Perez-Boerema, F., Pitto, L., Garcia Aznar, J. M., et al. (2020). A multi-scale modelling framework combining musculoskeletal rigid-body simulations with adaptive finite element analyses, to evaluate the impact of femoral geometry on hip joint contact forces and femoral bone growth. PLoS ONE 15, e0235966. doi:10.1371/journal.pone.0235966

PubMed Abstract | CrossRef Full Text | Google Scholar

Kainz, H., Killen, B. A., Van Campenhout, A., Desloovere, K., Garcia Aznar, J. M., Shefelbine, S., et al. (2021a). ESB Clinical Biomechanics Award 2020: Pelvis and hip movement strategies discriminate typical and pathological femoral growth – insights gained from a multi-scale mechanobiological modelling framework. Clin. Biomech. 87, 105405. doi:10.1016/j.clinbiomech.2021.105405

CrossRef Full Text | Google Scholar

Kainz, H., Wesseling, M., and Jonkers, I. (2021b). Generic scaled versus subject-specific models for the calculation of musculoskeletal loading in cerebral palsy gait: Effect of personalized musculoskeletal geometry outweighs the effect of personalized neural control. Clin. Biomech. 87, 105402. doi:10.1016/j.clinbiomech.2021.105402

CrossRef Full Text | Google Scholar

Kandzierski, G., Matuszewski, Ł., and Wójcik, A. (2012). Shape of growth plate of proximal femur in children and its significance in the aetiology of slipped capital femoral epiphysis. Int. Orthop. (SICOT) 36, 2513–2520. doi:10.1007/s00264-012-1699-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Kurz, M. J., Arpin, D. J., and Corr, B. (2012). Differences in the dynamic gait stability of children with cerebral palsy and typically developing children. Gait Posture 36, 600–604. doi:10.1016/j.gaitpost.2012.05.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Linde, F., Hvid, I., and Jensen, N. Chr. (1985). Material properties of cancellous bone in repetitive axial loading. Eng. Med. 14, 173–177. doi:10.1243/EMED_JOUR_1985_014_042_02

PubMed Abstract | CrossRef Full Text | Google Scholar

Maas, S. A., Ellis, B. J., Ateshian, G. A., and Weiss, J. A. (2012). FEBio: Finite elements for Biomechanics. J. Biomech. Eng. 134, 011005. doi:10.1115/1.4005694

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackay, J., Thomason, P., Sangeux, M., Passmore, E., Francis, K., and Graham, H. K. (2021). The impact of symptomatic femoral neck anteversion and tibial torsion on gait, function and participation in children and adolescents. Gait Posture 86, 144–149. doi:10.1016/j.gaitpost.2021.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Mellon, S. J., and Tanner, K. E. (2012). Bone and its adaptation to mechanical loading: A review. Int. Mater. Rev. 57, 235–255. doi:10.1179/1743280412Y.0000000008

CrossRef Full Text | Google Scholar

Mirtz, T. A., Chandler, J. P., and Eyers, C. M. (2011). The effects of physical activity on the epiphyseal growth plates: A review of the literature on normal physiology and clinical implications. J. Clin. Med. Res. 3, 1–7. doi:10.4021/jocmr477w

PubMed Abstract | CrossRef Full Text | Google Scholar

Modenese, L., and Renault, J.-B. (2021). Automatic generation of personalised skeletal models of the lower limb from three-dimensional bone geometries. J. Biomech. 116, 110186. doi:10.1016/j.jbiomech.2020.110186

PubMed Abstract | CrossRef Full Text | Google Scholar

Modenese, L., Barzan, M., and Carty, C. P. (2021). Dependency of lower limb joint reaction forces on femoral version. Gait Posture 88, 318–321. doi:10.1016/j.gaitpost.2021.06.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrell, D. S., Pearson, J. M., and Sauser, D. D. (2002). Progressive bone and joint abnormalities of the spine and lower extremities in cerebral palsy. RadioGraphics 22, 257–268. doi:10.1148/radiographics.22.2.g02mr19257

PubMed Abstract | CrossRef Full Text | Google Scholar

Peeters, N., Hanssen, B., Bar-On, L., De Groote, F., De Beukelaer, N., Coremans, M., et al. (2023). Associations between muscle morphology and spasticity in children with spastic cerebral palsy. Eur. J. Paediatr. Neurol. 44, 1–8. doi:10.1016/j.ejpn.2023.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Pizzolato, C., Lloyd, D. G., Sartori, M., Ceseracciu, E., Besier, T. F., Fregly, B. J., et al. (2015). Ceinms: A toolbox to investigate the influence of different neural control solutions on the prediction of muscle excitation and joint moments during dynamic motor tasks. J. Biomech. 48, 3929–3936. doi:10.1016/j.jbiomech.2015.09.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Rauch, F. (2005). Bone growth in length and width: The yin and yang of bone stability. J. Musculoskelet. Neuronal Interact. 5, 194–201.

PubMed Abstract | Google Scholar

Rethlefsen, S. A., Blumstein, G., Kay, R. M., Dorey, F., and Wren, T. A. L. (2017). Prevalence of specific gait abnormalities in children with cerebral palsy revisited: Influence of age, prior surgery, and gross motor function classification system level. Dev. Med. Child. Neurol. 59, 79–88. doi:10.1111/dmcn.13205

PubMed Abstract | CrossRef Full Text | Google Scholar

Rho, J. Y., Ashman, R. B., and Turner, C. H. (1993). Young’s modulus of trabecular and cortical bone material: Ultrasonic and microtensile measurements. J. Biomech. 26, 111–119. doi:10.1016/0021-9290(93)90042-D

PubMed Abstract | CrossRef Full Text | Google Scholar

Robin, J., Graham, H. K., Selber, P., Dobson, F., Smith, K., and Baker, R. (2008). Proximal femoral geometry in cerebral palsy: A population-based cross-sectional study. J. Bone Jt. Surg. Br. Vol. 90-B, 1372–1379. doi:10.1302/0301-620X.90B10.20733

PubMed Abstract | CrossRef Full Text | Google Scholar

Schless, S.-H., Hanssen, B., Cenni, F., Bar-On, L., Aertbeliën, E., Molenaers, G., et al. (2018). Estimating medial gastrocnemius muscle volume in children with spastic cerebral palsy: A cross-sectional investigation. Dev. Med. Child. Neurol. 60, 81–87. doi:10.1111/dmcn.13597

PubMed Abstract | CrossRef Full Text | Google Scholar

Shefelbine, S. J., and Carter, D. R. (2004). Mechanobiological predictions of femoral anteversion in cerebral palsy. Ann. Biomed. Eng. 32, 297–305. doi:10.1023/B:ABME.0000012750.73170.ba

PubMed Abstract | CrossRef Full Text | Google Scholar

Steele, K. M., DeMers, M. S., Schwartz, M. H., and Delp, S. L. (2012). Compressive tibiofemoral force during crouch gait. Gait Posture 35, 556–560. doi:10.1016/j.gaitpost.2011.11.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Stevens, S. S., Beaupré, G. S., and Carter, D. R. (1999). Computer model of endochondral growth and ossification in long bones: Biological and mechanobiological influences. J. Orthop. Res. 17, 646–653. doi:10.1002/jor.1100170505

PubMed Abstract | CrossRef Full Text | Google Scholar

Szuper, K., Schlégl, Á. T., Leidecker, E., Vermes, C., Somoskeöy, S., and Than, P. (2015). Three-dimensional quantitative analysis of the proximal femur and the pelvis in children and adolescents using an upright biplanar slot-scanning X-ray system. Pediatr. Radiol. 45, 411–421. doi:10.1007/s00247-014-3146-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhlrich, S. D., Jackson, R. W., Seth, A., Kolesar, J. A., and Delp, S. L. (2022). Muscle coordination retraining inspired by musculoskeletal simulations reduces knee contact force. Sci. Rep. 12, 9842. doi:10.1038/s41598-022-13386-9

PubMed Abstract | CrossRef Full Text | Google Scholar

van Arkel, R. J., Modenese, L., Phillips, A. T. M., and Jeffers, J. R. T. (2013). Hip abduction can prevent posterior edge loading of hip replacements. J. Orthop. Res. 31, 1172–1179. doi:10.1002/jor.22364

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Krogt, M. M., Bar-On, L., Kindt, T., Desloovere, K., and Harlaar, J. (2016). Neuro-musculoskeletal simulation of instrumented contracture and spasticity assessment in children with cerebral palsy. J. NeuroEng. Rehabil. 13, 64. doi:10.1186/s12984-016-0170-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Veerkamp, K., Schallig, W., Harlaar, J., Pizzolato, C., Carty, C. P., Lloyd, D. G., et al. (2019). The effects of electromyography-assisted modelling in estimating musculotendon forces during gait in children with cerebral palsy. J. Biomech. 92, 45–53. doi:10.1016/j.jbiomech.2019.05.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Veerkamp, K., Kainz, H., Killen, B. A., Jónasdóttir, H., and van der Krogt, M. M. (2021). Torsion Tool: An automated tool for personalising femoral and tibial geometries in OpenSim musculoskeletal models. J. Biomech. 125, 110589. doi:10.1016/j.jbiomech.2021.110589

PubMed Abstract | CrossRef Full Text | Google Scholar

Yadav, P., Shefelbine, S. J., and Gutierrez-Farewik, E. M. (2016). Effect of growth plate geometry and growth direction on prediction of proximal femoral morphology. J. Biomech. 49, 1613–1619. doi:10.1016/j.jbiomech.2016.03.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Yadav, P., Shefelbine, S. J., Pontén, E., and Gutierrez-Farewik, E. M. (2017). Influence of muscle groups’ activation on proximal femoral growth tendency. Biomech. Model Mechanobiol. 16, 1869–1883. doi:10.1007/s10237-017-0925-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Yadav, P., Fernández, M. P., and Gutierrez-Farewik, E. M. (2021). Influence of loading direction due to physical activity on proximal femoral growth tendency. Med. Eng. Phys. 90, 83–91. doi:10.1016/j.medengphy.2021.02.008

PubMed Abstract | CrossRef Full Text | Google Scholar