1 Introduction

Growth plate mechanobiology is a growing field of research that examines the mechanisms that regulate bone growth and development.^1–3 Located at the ends of long bones, growth plates drive longitudinal growth and are influenced by mechanical and biological factors.^4,5 Growth plates change from flat and smooth at birth to wavy with age due to mechanical stresses.^6,7,8 This shape change may prevent failure by shear stresses, as in slipped capital femoral epiphysis (SCFE).⁹ The thickness remains constant during childhood due to a balance between chondrocyte proliferation and apoptosis.¹⁰ Puberty disrupts this balance, initiating growth plate fusion and ossification, typically ending longitudinal growth around ages 17–19,¹¹ influenced by hormonal signals. Unlike humans, animals like mice and rats do not fully close their growth plates during maturation, yet they serve as a model of the influence of mechanical factors in bone development.^12,13

Diseases affecting the growth plate, such as Achondroplasia, Gaucher, and hip dysplasia, have significant prevalence, impacting patients' quality of life.^14–16 Salter-Harris fractures, constituting 30 % of pediatric fractures, can lead to growth disturbances if untreated.¹⁷ Animal models have been used to study physeal injuries, tissue engineering applications, and molecular mechanisms (Chen et al., 2009; Sundararaj et al., 2015.¹⁸; As an example, Wang et al., 2021 reviews different regeneration alternatives where different animal models are summarized and the importance of tissue engineering is highlighted as a main research area in the incoming years. In contrast, Mathematical modeling, enhanced by accessible numerical data and computational power, relates mechanical states of the growth plate to its development.¹⁹

Classic studies, like,²⁰ connected cartilage formation to stress types, using the osteogenic index (OI), extensively in bone development modeling, where a positive index means that there will likely be growth, and a negative index means inhibition or growth.Guevara-Morales et al. ²¹ analyzed stress impacts on growth plate shape, while²² studied tension-band plates' effects on knee growth plates using personalized FEM models, Yadav et al.²³ examined growth plate geometry's influence on femoral head morphology using a strain tensor based on the OI. These studies enhance understanding of the biomechanical environment role in skeletal development, suggesting biomechanical principles integration could lead to new treatments for growth disorders.

Advances in imaging and computational modeling reveal biomechanical forces and microstructural changes in the growth plate during growth and disease progression.^24,25 The growth plate's shape affects its load-sustaining capacity, especially under high shear stresses.^9,26 No current model predicts the growth plate shape change and its adjacent trabecular structures. This work introduces a new endochondral growth model focusing on trabecular patterns in growth plate development, aiding in understanding growth plate diseases and developing treatments targeting mechanical factors at different developmental stages. This model uses a bone remodeling approach based on strain energy density.^27,28 For cartilaginous zones, the evolution law is based on the OI,²⁰ determining strain rate related to proliferation and hypertrophy in the epiphyseal plate. This is the first use of a bone remodeling scheme to investigate localized load transmission's effect on growth plate shape and remodeling.

2 Methods

In this section, we present a modified bone remodeling algorithm, which incorporates a strain rate tensor dependent on the OI, which accounts for the net growth rate in the growth plate. The model considers a linear isotropic behavior for both bone and cartilage, this assumption is sustained in the fact that small displacements are considered. As a first step, a benchmark geometry will be addressed, then the model we propose will be applied to a pathologic case and healthy growth.

2.1

2.1 Model description

The remodeling algorithm shown in Fig. 1, was modified from a previous work where the relation between the shape of the growth plate and the trabecular patterns was analyzed via a static model that only considered bone remodeling,²⁹ whereas now, we consider the growth due to hypertrophy and proliferation in the growth plate through imposing a strain rate that depends on the OI in the lower elements of the growth plate. This new model also uses the strain energy density as the main stimulus that controls the bone remodeling process, as proposed by.^28,30–32 The elastic problem is defined by Eq. (1), which establishes the equilibrium between the divergence of the stress field σ and the body forces (b). This equation is coupled with the evolution law presented in Eq. (2), where the dimensionless density (λ) depends on the strain energy per unit volume (U) at the finite element, along with a reference strain energy (Uref), which serves as the threshold for bone formation or resorption. Additionally, the parameters (k1) and (n) represent the constant of remodeling speed and an experimental exponent, respectively. The remodeling and growth process are iterative, and at each step the elastic modulus (E) is updated following the power law shown in Eq. (3) until a convergence criterion is achieved.Eq. 1∇Tσ+b=0Eq. 2dλdt=k1[λn−1UstrUref−1]Eq. 3E(ρ)=E0λn

Close

Fig. 1 Modified bone remodeling algorithm incorporating growth plate dynamics. The algorithm predicts the formation of ossification bridges as a result of mechanical stimuli and bone adaptation with changes in the growth plate topology.

Export to PPT

The strain rate tensor ε˙ we proposed, follows a linear relationship with the OI which is added to account for the growth due to chondrocyte hypertrophy and proliferation in the n direction, (along the bone longitudinal axis) as seen in Eq. (4). The osteogenic index as proposed by³³ is a linear combination of the peak octahedral shear stress and the hydrostatic stress (Eq. (5)). The linear relationship between the strain tensor rate ε˙ and OI, has been previously explored by^23,34 which is useful to account for longitudinal growth. In particular, our model is based on a fitting of the growth rate in the experimental works of Villemure et al.³ who stated the contribution of growth rate in the hypertrophic and proliferative zones in a given time under normal loading conditions. Nonetheless, it is important to note that the daily elongation may vary for different bones. The parameters used throughout the simulations for bone remodeling in this study are seen in Fig. 1, which are taken from the works of.³⁰

The convergence criterion for the model was set to 150 days, with a timestep of 0.1 days, since density shows no visible changes after this point in the domain topology, similar to the work of Buenzli et al.³⁵ First, a quasi-stationary state is achieved for bone remodeling in the trabecular portion of the domain (until no appreciative changes can be seen, usually at 100 days), then at this point, the lower part of the growth plate in our algorithm is set to undergo remodeling, in the ossification zone, giving place to a bony connection in the growth plate, giving place to ossification bridges, seen at 150 days, numerically.Eq. 4ε˙=ε˙prolif+ε˙Hypertro=k2OIn⊗nEq. 5OI=τoct+kOIσhyd

2.2

2.2 Numerical implementation

The different simulations were solved using a user element subroutine (UEL, ABAQUS- 2017). The bone density evolution equation (Eq. (2)) is solved using an integration scheme of Runge Kutta 4th order. All simulations were done using a plane stress formulation, and linear triangular and quadrilateral elements were used, similar to the benchmark cases used by.^28,30 The mesh employed a characteristic length of 1 mm for each element. Since this type of topological problem is highly dependent on mesh size, a convergence analysis was performed for each mesh configuration, ensuring a change no higher than 5 % in stress. As a first step, a simple benchmark test used in other remodeling studies was proposed to represent a domain with both trabecular bone and cartilage properties from the growth plate. The bone properties were adopted from³⁶ and cartilage properties were averaged from those reported by.³⁷ In these benchmark tests, initial conditions correspond to two different densities for trabecular bone and cartilage, and an initial topology for the growth plate. Boundary conditions consist of a non-uniform load with supports, as seen in Fig. 2 (a). Different tests were performed using different topologies, i.e. a flat shape for the growth plate with different inclinations, concave form (which is typical in advanced stages of development in the femur epiphyseal plate) or a sine wave pattern describing the shape of the growth plate.³⁸ Each configuration has a distinct set of trabecular formations and corresponding growth plate responses. In the annexes, various shapes and their resulting topologies are presented, along with a surface response to evaluate effect of the growth plate shape on bone density. This analysis reveals that the primary factors influencing bone density along the growth plate are thickness (p = 0.03) and the number of oscillations, which shows marginal significance (p = 0.05).

Close

Fig. 2 a) Different trabecular adaptations for different topologies of the growth plate, in the center initial flat shape, on the right, concave shape (growth plate radius of 0.3m).

Export to PPT

These benchmark tests were used to calibrate the model and to perform convergence analysis. Once the model was calibrated, a hip dysplasia case, and the growth of a proximal tibia were studied to appreciate the effect of abnormal loading during development and its dependency on load and trabecular and growth plate cartilage topology. The boundary conditions were adapted from Peyroteo et al.³⁹ to model the extreme range of abduction and adduction in the femur, and finally,Torcasio et al.⁴⁰ for the compression in the tibial plateau.

3 Results

The results demonstrate how changes in growth plate shape can emerge as a consequence of bone remodeling in its vicinity, while trabecular structures adapt in response to mechanical stimuli. This process is illustrated in Fig. 2(a), which presents the implementation of the algorithm on a benchmark test with a predefined growth plate shape. Over time, distinct morphological changes occur: after 100 days, the algorithm allows the remodeling of the lower elements in the growth plate's midsection. By 150 days, structures resembling ossification bridges begin to form at the periphery, accompanied by the emergence of mammillary process-like formations in the surrounding bone, as shown in Figs. 2(a) and 4. Additional simulations displaying different growth plate topologies and their corresponding mammillary processes are included in the annexes.

A regression analysis was conducted to quantify the relationship between growth plate morphology and bone density (p = 0.183). The analysis identified thickness and amplitude as the most influential factors affecting bone density near the growth plate. The response surface, illustrated in Fig. 3, suggests a high degree of nonlinearity, with strong interaction effects between thickness and amplitude (statistical details are provided in the supplementary material). A Spearman correlation analysis further confirmed these trends, revealing the strongest positive correlation for thickness (r = 0.38), followed by amplitude (r = 0.36). In contrast, frequency exhibited a negative correlation (r = −0.33), indicating that higher frequencies may lead to lower bone density.

Close

Fig. 3 a) Surface response of Bone Density vs. Amplitude & Thickness. Frequency is represented by the diameter of each data point.

Export to PPT

Close

Fig. 4 Benchmark test showing evolution of the growth plate and the trabecular patterns in its vicinity, growth plate bridging and similar structures to mammillary processes are seen.

Export to PPT

Building on the hypothesis that trabecular patterns influence growth plate morphology, we investigated the role of the osteogenic index (OI) as a predictor of ossification bridge formation in a clinical case of hip dysplasia (HD). This condition is characterized by altered mechanical loading due to reduced lateral coverage of the femoral head, which we simulated using the boundary conditions shown in Fig. 5(c). Specifically, we applied a non-uniform radial load, typical of HD, with a joint reaction force of 1158 N (Beaupre and Orr 1990) at an inclination of 34.5°, corresponding to the pathological neck shaft angle depicted in Fig. 4(a).

Close

Fig. 5 a) Hip dysplasia in a skeletally mature individual, modified from²⁵ b) Trabecular groups and enlarged femoral head for flat topology (left) and wavy pattern in the growth plate (right). c) Boundary conditions and OI in developmental dysplasia matching trabecular patterns.

Export to PPT

To assess the influence of growth plate geometry on load transmission and subsequent trabecular adaptation, we tested two distinct topologies: a flat growth plate and a wavy-patterned growth plate (Fig. 5 (b)). In both cases, the altered loading conditions induced significant changes in trabecular architecture, particularly in the compression groups, leading to a reduction in the neck shaft angle and elongation of the femoral head (Fig. 5 (c)). Notably, these modifications resulted in non-uniform growth modulation, with localized zones of growth arrest contributing to deviations in shaft angle and neck length. Despite these variations, the final trabecular topology remained consistent across both growth plate configurations.

Notably, the formation of ossification bridges correlated with the osteogenic index, reinforcing findings from the benchmark cases. This suggests that localized mechanical stimuli play a crucial role in pathological bone adaptation. Understanding these localized growth modulations could be relevant for developing treatment strategies for hip dysplasia, potentially allowing for targeted interventions based on mechanical load distribution.

Fig. 6 shows the tibial plateau of a C57BL/6 mouse at two and three months, highlighting the formation of ossification bridges and the associated changes in trabecular patterns. These trabecular structures align with the bridges, suggesting a mechanical link between growth plate morphology and trabecular adaptation. Additionally, variations in growth plate thickness and amplitude can be observed over time.

Close

Fig. 6 (a) Tibial plateau at two months (top), showing the absence of ossification bridges, compared to 21 months (bottom), where ossification bridges are present. (b) Algorithm applied to a reconstructed 2D geometry, illustrating longitudinal growth (see green line), which totals 396 μm.

Export to PPT

Simulations of the tibial plateau illustrate that ossification bridges primarily form at the periphery and at the inflection points of the growth plate curvature. This trend is consistent with early-stage trabecular organization, where trabeculae are more uniformly distributed near the growth plate, compared to later stages when the topology becomes more complex and oscillations start to appear along with a diminished thickness over time (Fig. 5(a)(b)). This pattern likely results from increased growth plate width, which helps distribute stresses more evenly (see Fig. 2,Fig. 6). The statistical analysis from the benchmark case supports this observation, as thickness emerged as the strongest coefficient influencing trabecular adaptation (Surface fit, and Pearson coefficients shown in the Annexes).

4 Discussion

This study presents a computational model that utilizes the osteogenic index (OI) to predict growth plate deformation over time, considering both the dynamic behavior of trabecular bone and growth plate cartilage. The model simulates the biomechanical environment, considering a linear elastic model, incorporating bone remodeling and the development of trabecular patterns and endochondral growth in response to mechanical loads, it adds on previous work where only remodeling is considered.²⁹ Our numerical results indicate that growth plate shape significantly influences bone density in its vicinity, alters trabecular architecture on both distal and proximal sides, and may serve as an adaptive mechanism to better withstand shear stress, with bone density adaptations, similar to mammillary processes. These findings align with previous cadaveric studies⁹ and support the broader biomechanical principle that bones respond to mechanical stimuli through structural adaptation, as described by Wolff's law (Cowin, and Telega, 2003;⁴¹). Our approach extends this principle by explicitly modeling the interaction between growth plate cartilage and trabecular bone.

Initially, the model's ability to generate trabecular patterns and its interaction with the growth plate were evaluated. This was followed by an assessment of overall longitudinal growth, identifying distinct zones of growth modulation, nevertheless, here we assume that growth is either promoted or inhibited according to the sign of the OI, a future work could define thresholds where lazy zones or modulated growth can occur, depending on the magnitude and frequency of the mechanical stimuli,⁴² furthermore other studies have defined a threshold effect for chondral modelling.⁴³ Lastly, the formation of ossification bridges was analyzed. These aspects were qualitatively validated by comparing different growth plate morphologies and examining bone formation in regions of high shear stress (Fig. 6(a)–(b)). To extend this analysis, a quantitative framework was developed to assess the influence of growth plate shape on bone formation and growth. Benchmark tests were conducted to examine bone formation along the growth plate and growth rates, which followed a linear relationship with the osteogenic index (OI), a widely used parameter in cartilage mechanobiology.²⁰

These findings suggest a two-way biomechanical interaction, where growth plate morphology impacts trabecular architecture, and in turn, trabecular adaptations influence growth plate shape evolution.The way that bone may help maintain mechanical stability has been addressed before by,²⁶ who analyzed the development of mammillary processes in the porcine femoral head, proposing that their radial pattern enhances biomechanical stability and may play a role in preventing slipped capital femoral epiphysis (SCFE). Notably, our algorithm updates the strain rate tensor at each iteration after recalculating bone density (λ), ensuring that mechanical adaptations are dynamically incorporated. Additional simulations, including baseline profiles of τoct, σhyd, and OI, are available in the supplementary material. As expected, high-density zones act as stress concentrators, which further enhance local density, and changes in curvature coincide with peaks in shear stress values, reinforcing the mechanical link between these factors.

This model has significant clinical applications, particularly for guided growth procedures and the correction of bone deformities. It could also be used to assess pathologies such as cam deformity, which are associated with mechanical forces.⁴⁴ Hemiepiphysiodesis, a surgical technique that temporarily arrests growth in a specific region of the growth plate to correct angular deformities, is an example of growth modulation due to mechanical stress playing a crucial role. Ding et al. ⁴⁵ investigated the histopathological effects of hemiepiphysiodesis in a mini pig model, demonstrating that asymmetric pressure distribution led to differential growth plate thickness changes. Their study highlighted that stress concentration in specific regions induced apoptosis and altered the expression of key molecular markers such as Caspase-3, Caspase-9, P65, and CHOP, which are involved in chondrocyte apoptosis and endochondral ossification. Notably, the removal of the implant resulted in a thickening of the growth plate in previously compressed regions, contributing to the recurrence of malformations. To investigate further cases such as this one, our model provides a framework to investigate the biomechanical factors influencing growth arrest, offering valuable insights into optimizing surgical interventions and predicting post-surgical outcomes.

Beyond growth modulation, conditions such as congenital pseudarthrosis, a disorder associated with abnormal trabecular patterns in the growth plate, may also benefit from this mechanistic approach. Understanding the formation and role of trabecular patterns could contribute to better therapeutic strategies for disorders including SCFE, Legg-Calvé Perthes disease and Hip Hysplasia. Simulations of hip dysplasia and trabecular adaptation reveal a direct relationship between growth plate morphology, trabecular architecture, and overall growth rate that may lead to deformities in the femoral head. Conditions where shear stress is a key factor, such as slipped capital femoral head and Sever's disease, could also benefit from similar modeling strategies. Moreover, documented variations in growth plate waviness in these conditions further support the need for biomechanical models that account for stress-dependent changes in growth plate topology.⁸

Our μCT analysis revealed a strong correlation between high shear stress zones, often located at growth plate inflection points and the initiation of ossification bridges. A finite element study by Castro et al.³⁸ further demonstrated that in a wavy-patterned growth plate, strain decreases as irregularities increase. This expectedly amplifies stress, making pronounced inflection points susceptible to tissue failure.

Ossification bridges, which emerge at the onset of ossification, connect trabecular groups above and below the epiphyseal growth plate, particularly in regions of elevated shear stress. This pattern is evident in Fig. 7(c), where ossification bridges appear at inflection points, coinciding with peak shear stress and linking trabecular structures, similar to observations in Fig. 7 (a) and 7 (b). Additionally, Fig. 7 (d) illustrates ossification bridges in a 20-month-old C57BL6 mouse, reinforcing this trend. Notably, these bridges are most prominent in the peripheral growth plate, where shear stress is highest, aligning with findings from Staines et al. (2018) on the tibial plateau.

Close

Fig. 7 a) Change in width of the growth plate of a C57BL/6 mice and formation of ossification bridges (shown in yellow). b) Rib from a 12- week-old female infant⁴⁶ showing high density zones. c)Trabecular formations for a wavy growth plate, bridge formation at 150 days in yellow. (d) Ossification points observed at three months, predominantly located in regions of high shear stress.

Export to PPT

Simulations of hip dysplasia and trabecular adaptation reveal a direct relationship between growth plate morphology, trabecular architecture, and overall growth rate uniformity, suggesting potential applications in growth prediction models. A 2D computational framework could be extended to incorporate the stiffness and placement of staples or screws, allowing for a more precise prediction of growth modulation, in our model we proposed a growth modulation based on the osteogenic index. Conditions where shear stress plays a significant role, including slipped capital femoral head and Sever's disease, could also benefit from similar modeling approaches. Furthermore, documented changes in growth plate waviness in these conditions reinforce the need for biomechanical models that account for stress-dependent variations in growth plate topology.

In summary, this study presents a novel computational approach to understanding the biomechanical interplay between growth plate deformation and trabecular bone remodeling. By integrating both qualitative observations and quantitative modeling, we demonstrate that growth plate shape influences trabecular adaptation, particularly in response to shear stress, and that ossification bridges form at high-stress regions, affecting skeletal development, furthermore we identified, numerically, the thickness and the amplitude of the growth plate as main parameters influencing bone density. The integration of mechanobiological factors into growth plate modeling could enhance treatment strategies for growth-related disorders, improving the precision of intervention planning in pediatric orthopedics. Future work should focus on refining this model by incorporating anisotropic material properties, patient-specific geometries, and 3D growth plate dynamics, further bridging the gap between computational predictions and clinical applications.

5 Conclusions

This study presents a model that integrates bone remodeling and endochondral growth to examine the evolution of growth plate shape during development. A key aspect of this model is its ability to account for trabecular patterns in the surrounding bone, which interact with the growth plate and influence its morphology over time. The generated topologies were validated by comparison with clinical data and μCT scans from mice. Additional simulations, along with a statistical analysis of the effect of growth plate shape on bone density, are provided in the supplementary material, further supporting these findings. This validation establishes a correlation between regions prone to ossification, bridge formation and alterations in growth plate topology, offering insights into the development of mammillary processes, structural features critical for bone stability.

The model relies on various assumptions, particularly regarding material properties. Both bone and cartilage were modeled as linear isotropic elastic materials, despite cartilage being known for being a poro-viscoelastic-fibribilar⁴⁷ material and bone having an anisotropic behavior. This simplification is justified by the small strain conditions considered, yielding results consistent with existing medical literature on bone and epiphyseal cartilage morphology, specifically we were able to obtain different anatomical features such as trabecular groups, localized change in the growth plate topology, and overall bone growth rates. Furthermore, the strain tensor rates produced minimal nodal displacements relative to element size, emphasizing the critical role of the initial growth plate topology in determining morphological outcomes. Future improvements could involve remeshing techniques within the user subroutine to better capture morphological changes over time.

A potential application of this model is in guided growth therapies, such as hemiepiphysiodesis or transphyseal screw-based growth modulation. By quantifying growth modulation across the entire growth plate and tracking geometric changes, the model could support the optimization of therapeutic strategies. Additionally, studying the adaptive response of the growth plate to mechanical stimuli has implications for diagnosing and treating growth deformities, including angular deformities and conditions like slipped capital femoral epiphysis. The ability to dynamically analyze trabecular patterns may also enable early disease detection, facilitating timely interventions.

Beyond these applications, this research provides a computational framework that links bone remodeling and endochondral growth, offering a mechanistic perspective on growth plate adaptation. The model has the potential to support clinical decision-making by simulating how mechanical and biological factors shape growth plate morphology. Future refinements, such as incorporating anisotropic material properties and patient-specific geometries, could further improve its predictive accuracy and clinical relevance.

Author contribution

The study was collaboratively designed by all authors as part of a cotutelle Ph.D. program between the Universidad Nacional de Colombia and the Université de Technologie de Compiègne. The first author, Diego Alfredo Quexada Rodríguez, drafted the manuscript, with the assistance of Professors Olfa Trabelsi (maître de conference)from Université de Technologie de Compiègne and Diego Garzón (full professor), who provided valuable corrections. Professor Marie-Christine Ho-Ba-Tho (full professor) from Université de Technologie de Compiègne played a crucial role in conceiving the study and providing insightful guidance. Software used in this study consists of ABAQUS (2018), Matlab (2023). All funding was acquired thanks to Université de Technologie de Compiègne and Universidad Nacional de Colombia as part of the aforementioned cotutelle program. No conflict of interest exists.

Ethical statement

All experimental procedures involving human subjects, animals, or sensitive biological materials have been conducted in accordance with ethical guidelines established by relevant institutional review boards (IRBs) or ethics committees. Necessary approvals and informed consent statements have been obtained and documented.

Patient consent

The present study did not involve any prospective data collection from human participants or identifiable personal information.

Funding

This research was fully funded by the cotutelle PhD program established between Université de Technologie de Compiègne (France) and Universidad Nacional de Colombia (Colombia). The funding supported the collaborative research, mobility, and resources necessary for conducting this study as part of the doctoral training framework.