1 Introduction

Bones are mineralized connective tissues that provide structural support, protect organs, and enable movement.^1,2 At the organ level, bones consist of cartilaginous joints, growth plates, marrow spaces, and mineralized cortical and trabecular structures,³ at the tissue level, they are composed of a cellular matrix that includes osteoblasts, osteocytes, and osteoclasts, which are essential for skeletal growth, modeling, and remodeling.^2,4 Osteoblasts form new bone tissue by synthesizing and secreting the bone matrix, contributing to mineralization. Osteocytes, derived from osteoblasts, maintain the bone matrix and regulate mineral content. In contrast, osteoclasts are multinucleated cells responsible for bone resorption, playing a crucial role in bone remodeling and healing.²

While osteocytes and osteoclasts maintain and remodel existing bone, the formation of new bone relies on other cell types and a process called endochondral ossification, which is the process responsible for the formation of most long bones, such as the femur and tibia.⁵ This process begins begins with the condensation of mesenchymal cells, which initially differentiate into chondrocytes and produce a cartilage matrix, where new bone will develop.^6,7 The initial site of bone formation is the primary ossification center (POC), which emerges within the diaphysis or central shaft region of the long bone. Subsequently, additional sites of bone formation, referred to as secondary ossification centers (SOC), develop at the ends of the bones, or the epiphyses.⁸ The primary and secondary ossification centers are separated by the growth plate, which facilitates the longitudinal growth of the bone.^9,10 This growth plate consists of several zones: the resting zone with hyaline cartilage, the proliferation zone where chondrocytes undergo mitosis, and the hypertrophic zone where chondrocytes enlarge, promoting stiffness and calcification.^6,9 Vascular invasion then allows osteoblasts to form woven bone that remodels into lamellar bone, forming the trabecular structures essential for bone strength.^3,9

The rate and time of appearance of ossification centers is essential for evaluating skeletal maturity and recognizing growth-related conditions, such as limb length discrepancies and fractures involving the epiphyseal growth plates.^11,12 This knowledge provides valuable insights that support the diagnosis, management, and treatment of various developmental and traumatic musculoskeletal disorders. Time and rate of bone ossification can vary among different bones, and the process may also differ between genders.¹¹

The appearance of secondary ossification centers (SOCs) marks a crucial stage in bone development. Their formation is significantly influenced by the mechanical environment and plays a protective role for the growth plate against mechanical loading, as observed across diverse vertebrate taxa.¹³ This is particularly relevant in the context of synovial joints—articulations where two or more bones meet. These bones are covered by hyaline cartilage, which originates from the embryonic cartilage anlage. The primary function of these joints is to facilitate movement, aided by a lubricating film of synovial fluid in between the articular surfaces, and a fibrous capsule to contain the joint.^14,15 Research on quadrupedal mammals has explored the scaling of elbow dimensions in relation to body mass and locomotion, revealing insights into joint lubrication mechanisms.¹⁶

In terms classification, of synovial joints can be categorized based on their complexity into three types: simple, where the joint consists of two articular surfaces in contact; compound, with three or more articular surfaces; and complex, with a meniscus between the articular surfaces (see Fig. 1).¹⁵ Consequently, the epiphysis of certain bone can be subjected to multiple load points, each one corresponding to the structure in contact to it. For instance, in the elbow joint we have the distal humerus subjected to loading in two zones one exerted by the radius ans the other by the ulna (Fig. 1B). Furthermore, the range of movement and place of the applied load is defined by the shape of the joint, the larger the range of movement of certain bone, the more extended the load is on the surface of the epiphysis.

Close

Fig. 1 Different Types of Synovial Joints: (A) Simple joint with a single synovial cavity, (B) Compound joint with multiple articulating surfaces, and (C) Complex joint with an articular disk or meniscus. Bones can exhibit various geometries, including both convex and concave epiphyseal surfaces. Taken from Gray's anatomy¹⁵.

Export to PPT

Up to this point, different theories have been proposed to explain endochondral ossification progression. Some suggest mechanical load environment is a key factor on this process,^5,17–26 proposing that cyclic octahedral shear stresses stimulate ossification, while cyclic compressive dilatational stresses inhibit it (Fig. 2), potentially due to various mechanobiological mechanisms. For instance, cyclic shear stress has been shown to increase collagen-I formation²⁷ while simultaneously decreasing the production of proteins related to cartilage matrix, whereas cyclic hydrostatic pressure promotes their synthesis.²⁸ Furthermore, cyclic shear stress increases nitric oxide production, which is associated with chondrocyte apoptosis; whereas cyclic hydrostatic pressure can decrease its production.²⁸ Hydrostatic pressure creates pressure gradients that promotes fluid flow and nutrient delivery supporting chondrocyte viability, whereas shear stress inhibits these processes.^5,29 Additionally, aggrecan matrix production is stimulated by cyclic hydrostatic pressure,^28,29 which plays a crucial role in mediating chondrocyte-matrix interactions.³⁰ Therefore, by combining cyclic hydrostatic and octahedral shear stresses, a parameter can be created to account for the likelihood of cartilage at a particular point transforming into bone, as expressed by the ossification index (OI), detailed in Equation. (5).¹⁷

Close

Fig. 2 The hydrostatic stress and the octahedral shear stress³¹.

Export to PPT

Given the established importance of mechanical influences on bone development, and considering that the synovial joint significantly shapes the mechanical environment within a developing epiphysis, the nature of the synovial joint can directly affect the development of the epiphysis and its SOCs. Nevertheless, biochemical processes governing chondrocyte behavior are equally essential. In particular, chondrocyte hypertrophy is a key driver of endochondral ossification. Several theories have been proposed to explain the regulation of this process, emphasizing the interplay between parathyroid hormone-related peptide (PTHrP) and Indian hedgehog homolog (Ihh). PTHrP delays chondrocyte hypertrophy, while Ihh promotes it.³² These signaling molecules regulate each other, establishing a feedback loop in which Ihh stimulates PTHrP production, and PTHrP inhibits Ihh secretion.^32,33 Additionally, hypertrophic chondrocytes contribute to ossification by promoting matrix mineralization, vascularization, and the recruitment of bone precursor cells.^34–36 These processes are mediated by the secretion of signaling molecules such as vascular endothelial growth factor (VEGF).^37–40 To better understand this regulatory mechanism, mathematical models based on reaction-diffusion frameworks have been developed to simulate the PTHrP-Ihh feedback loop and its effects on chondrocyte hypertrophy and ossification, often in combination with mechanical models.^{5,14,20,21,41} Furthermore, the process is under tight genetic control, with key transcription factors such as Sox9 and Runx2 orchestrating chondrocyte differentiation and hypertrophy.^42,43 Systemic factors like thyroid hormone and growth hormone, among other endocrine influences, also regulate chondrocyte behavior on bone development.^42–45

Conversely, in terms of geometric representation of synovial joints, computational models have approached their morphology in different ways. Some models focus on two cartilage rudiments set in contact, with analysis of their mechanical and biochemical interactions used to predict the resulting joint shape,^21,23,24 or by focusing on a single cartilage rudiment and modeling the influence of the adjacent one.^{5,13,17,19,20,46,47} These models often rely on a fixed geometry and load conditions, presenting limited options of analysis, typically considering at most three geometry variations or three loading conditions. Even though some models incorporate changes in shape,^20,21,23,24 these emerge as results of the system's evolution rather than being defined as initial conditions. This limitation restricts the ability to systematically analyze a wide range of geometries and loading scenarios, which could be crucial for developing more biofidelic models of epiphyses. To address this gap, the model proposed in this work allows both the geometry and loading conditions of the rudiment to be defined a priori. This enables the analysis of various shapes and the influence of specific geometric and loading parameters by isolating their effects on the ossification pattern, including the number, shape, and location of SOCs. This study presents a simplified analysis, solely focusing on the morphological and mechanical influence on the ossification process, within a simplified two-dimensional, linear-elastic model. Similar approaches have also been proven effective in providing insights about endochondral ossification, yielding results consistent with clinical and experimental observations. Therefore, we do not consider more complex models in terms of constitutive form, biofidelity or dimensionality. Our focus is primarily on the parametric modeling approach, that differs from previous studies that often relied on a limited number of geometry and load cases. Additionally, this study lays a reference framework for future models to incorporate biochemical interactions and refine the mechanical approach by adding poroelasticity, bone remodeling or growth. Then, by modifying the model parameters, normal and abnormal geometric and loading conditions can be explored, representing pathological conditions and helping to identify underlying causes. Also, the model herein proposed offers a novel approach to explore epiphyseal shape across diverse animal species, correlating the appearance of secondary ossification centers with mechanical loading and rudiment shape.

2 Methodology

The proposed methodology begins with a CAD-based geometric model representing cartilage and bone regions using NURBS curves, which offer flexibility through control points. Seven geometrical parameters, such as angles and lengths, are left free, while parabolically distributed pressures are applied along the articular surface, each defined by two parameters dependent on geometry (location and span) and one parameter related to load intensity. An additional parameter sets the maximum load magnitude, with all other loads expressed as proportions of this value. A two-dimensional parametric space, defined by variables p and q, generates combinations of geometries and loading conditions, with geometry-driven load assignment simplifying the analysis. These configurations are solved using the Finite Element Method under plane strain assumptions and isotropic material properties, with bilinear quadrilateral meshing and convergence analysis ensuring reliability. Finally, the ossification index (OI) is computed to predict regions likely to form secondary ossification centers (SOCs), with results aligning with known biological observations. A graphical summary of the methodology is shown in Fig. 3.

Close

Fig. 3 The proposed methodology consists of several key steps for analyzing bone morphogenesis. These include defining the geometry using NURBS curves, analyzing the parameter space by sampling the free geometric parameters, assigning loads based on the geometry, discretizing the problem, and solving it using the Finite Element Method. Additionally, a convergence analysis is conducted to ensure the reliability of the results.

Export to PPT

2.1

2.1 Constitutive model

The physical model consists of two material phases: bone and cartilage. Bone, represented by the region Ωb, and cartilage, represented by the region Ωc, are modeled using a linear-elastic, isotropic, single phase approximations. The union of Ωb and Ωc, i.e., Ω = Ωc ∪ Ωb constitutes the entire domain Ω (see Fig. 4). The boundary conditions applied consist of Dirichlet constraints at the bottom surface and pressure load distributions at the top, which are defined in Section 2.4.

Close

Fig. 4 The constitutive model consists of two adjacent regions representing bone and cartilage materials. Boundary conditions are defined, including a Dirichlet condition on the bottom to keep this region fixed, and a Neumann condition where a distributed loads are applied.

Export to PPT

For a point x = [x, y]^T within Ω, static equilibrium must be satisfied. This requirement leads to the conservation of momentum equation. Given that body forces are neglected, the resulting equation is:(1)∇⋅σ=0,x∈Ω

The relationship between stress and strain follows Hooke's law for continuous media:(2)σ=Dϵ

where, in Voigt notation, the stress and strain tensors are represented as 3-dimensional vectors, and D acts as the constitutive matrix for the isotropic 2D plane strain case:(3)D=E(1−2ν)(1+ν)1−νν0ν1−ν00012−ν

strain ɛ is defined in terms of derivatives of displacements u = [u, v]^T in space.(4)ϵ=∂u∂x∂v∂y∂u∂y+∂v∂x

The material properties for bone (Ωb) and cartilage (Ωc) are defined based on values reported in the literature. For bone, a Poisson's ratio of 0.2 and a Young's modulus of 500 MPa are used, which are typical for newly mineralized bone.^13,17,48 For cartilage, properties are based on bovine samples, which are commonly used in research due to their structural similarity to human cartilage and easy availability. The shear modulus of bovine cartilage has been measured between 2 and 3 MPa.^31,49 A value of 2.04 MPa is selected within this range. Cartilage is nearly incompressible, so a Poisson's ratio close to 0.5 is appropriate; here, 0.47 is used to ensure numerical stability. The Young's modulus is then calculated using the relation E = 2G(1 + ν), resulting in E ≈ 6 MPa.^17,18,31

Consequently, Equation (1) and these boundary conditions define a system of two Partial Differential Equations (PDEs) with two unknowns, displacement u and v defined by the previous equation and boundary conditions can be solved to provide the displacement, strain, and stress fields. Additionally, the general ossification index for each point in the domain can be computed by calculating the individual OI values for each load in the loading history and taking the mean value.(5)OI=1n∑i=1nSi+kDi

where, in terms of the principal stresses (σi(i = 1, 2, 3)), the octahedral shear stress is:(6)S=(σ1−σ2)2+(σ2−σ3)2+(σ3−σ1)23

and hydrostatic stress is:(7)D=σ1+σ2+σ33

The coefficient k determines the relative contribution of shear and hydrostatic stresses. A value of 0.5 has been suggested in previous studies as suitable for modeling postnatal ossification, as values between 0.3 and 1.0 provide a realistic balance between stress components while avoiding predictions that are inconsistent with observed skeletal development.^5,17,18

2.2

2.2 Geometric model

The methodology herein proposed is strongly based on a tight control of the geometrical configuration. For the purposes of this work the selected geometry consists of Non-uniform rational B-splines (NURBS) curves, whose definition is provided in Section A.2. For the current case, the control polygons of the curves are defined by the geometric constraints available in the computer-aided design (CAD) model, such as dimensions, coincident points, symmetries, and tangencies (Fig. 5). Moreover, the geometry of the model is characterized by seven key dimensional parameters. Three parameters describe the head: radius_x (rx), radius_y (ry) and head_angle (αh). Two parameters define the ossification front: curve_angle (αf) and cartilage_thickness (tc). The remaining two parameters characterize the diaphysis: the head_heignt (hh) and the bone_width (wb). For further details, the CAD model has been made available online.

Close

Fig. 5 Parameters of the long bone geometry. The model is symmetric along the y axis and different parameters define the sizes of the epiphysis and diaphysis.

Export to PPT

As shown in Fig. 6, the geometry is enclosed by a set of external and middle curves that define the physical boundaries of the bone and cartilage domains. The central horizontal curve represents the ossification front, marking the transition from cartilage to bone. Internal construction lines were added to subdivide the cartilage region into quadrangular patches, enabling a more regular meshing scheme. The bottom curve represents the fixed boundary condition applied in simulations.

Close

Fig. 6 The external and middle curves define the physical boundaries of cartilage and bone, while internal construction curves divide the cartilage region into quadrilateral areas for mesh generation. The ossification front separates the bone and cartilage domains, and the fixed boundary is applied at the bottom. The articular surface is the superior.

Export to PPT

An advantage of using NURBS curves is that they provide knot points as essential reference markers for modeling contact surfaces. By utilizing the locations of these knot points and the distances between points on the curve, we can describe the load location, as they closely relate to curvature and, consequently, to biological surfaces. For instance, for concave surfaces, knots positioned at the external edge of the concavity—where the curve transitions direction—consistently mark the position of this phenomenon across different curves. Further details on Section 2.4.

2.3

2.3 Parametric analysis setup

The parametric analysis explores a two-dimensional space of geometries, defined by two independent variables, p and q (see Fig. 7). Each point in this space represents a distinct configuration. Although the model involves seven geometric parameters, each one is expressed as a function of p or q, based on selected minimum and maximum values and the desired number of samples. In our implementation, each parameter varies linearly with one of the parametric variables or remains constant.

Close

Fig. 7 Joint morphology and contact conditions.

Export to PPT

We analyze three scenarios, each one with its own parametric space of geometries. In scenario 1, parameters such as rx, ry, and αh, are linearly dependent on p, while tc is linearly dependent on q. In the other scenarios, rx depends on p, and ry depends on q. This parametric framework allows for the exploration of 25 possible variations for each scenario. Later, integer samples of p and q are selected, which determine corresponding values for the dependent parameters (Table 1). The first scenario shows a continuous transition between the concave and convex head shapes, allowing analysis of cases where the contact surface is planar and parameters such as the advancement of the ossification front are varied. The other two cases analyze concave and convex head shapes respectively, while varying variables related to the morphology of the head. The model is updated for each combination of p and q, generating a set of geometries. The ranges of the parametric variables were selected to reproduce epiphysis morphologies reported in the literature.^17,19

Table 1 Sample values for the different geometry cases. The curve angle (αf) was held constant at 15°, the bone width (wb) at 2.2 mm, and the head height (hh) at 3 mm on all cases. Head Angle (αh) and Cartilage Thickness (tc) were also held constant on cases 2 and 3.

	Scenario 1					Scenario 2					Scenario 3
	Convex to concave head					Concave head					Convex head
	a)	b)	c)	d)	e)	a)	b)	c)	d)	e)	a)	b)	c)	d)	e)
rx (mm)	1.4	1.4625	1.525	1.5875	1.65	1.1	1.325	1.55	1.775	2.0	1.1	1.325	1.55	1.775	2.0
rx (mm)	1.5	1.425	1.35	1.275	1.2	2.0	1.775	1.55	1.325	1.1	2.0	1.775	1.55	1.325	1.1
αh (°)	25.0	15.0	5.0	−5.0	−15	−15					15
tc (mm)	1.5	2.0	2.5	3.0	3.5	2.0					2.0

2.4

2.4 Loading conditions

Boundary conditions for a general case were declared on Section 2.1, however, for this specific case, a more precise definition is required. The fixed boundary, denoted a Γu, corresponds to the lower surface, while the top surface, Γσ, is where the pressure is applied. Two distinct loading history cases are considered: one in which the bone is subjected to load directly from another structure, such as in a simple or complex synovial joint. While in the other case, the load is applied by two adjacent bones, as observed in a compound synovial joint (see Fig. 8). The first case will be referred to as the Single Contact Case, while the second will be named the Double Contact Case.

Close

Fig. 8 Joint morphology and contact conditions in various synovial joints. The top row illustrates examples of different synovial joints and the bones they connect, while the bottom row shows combinations of articulating surface geometries and associated loading conditions.

Export to PPT

In terms of loading, each point on the articular surface has a defined pressure value following a parabolic load distribution (Fig. 9). The distribution is characterized by three parameters: the vertex coordinates of the parabola (hi, ki) and a term relative to the area over which the pressure is applied (ri). The equation of a single load is:(8)P=(t−hi)24pi+kiwherepi=−ri24ki

Close

Fig. 9 Schematic representation of an individual load. Consider a coordinate system that extends along the articular surface. Its origin is located at the center, with the positive tangential coordinate (t) extending to the right and the normal coordinate (n) to this line at each point. Each individual load is defined as an inverted parabola, characterized by the maximum point (ki), horizontal vertex position (hi), and load radius (rl), which represents the total span of the load in the tangential direction.

Export to PPT

Knots are characteristic points on a NURBS curve (see section A.2 for the mathematical formulation). By construction, the curve along which the loads are applied has an odd number of knots, with the middle knot located on the y-axis (highlighted by the red circle in Fig. 9). The distance between adjacent knots enables the definition of three characteristic lengths used for load placement (illustrated on Fig. 9): Lmin, corresponding to the distance between immediately adjacent knots, Lmid the distance between the next pair, and Lmax, defined by the subsequent pair beyond.

The total force applied (Fi) can be calculated by integrating the pressure over the load span and multiplying by the thickness of the model on z-coordinate (T). Due to symmetry, the total force is given by:(9)Fi=2kiT∫hihi+ri1−t−hiri2dt

Performing the integration and simplifying the expression yields:(10)Fi=43kiriT

From this relationship, we derive an expression for ki as:(11)ki=3Fi4riT

Previous studies have assumed that joint contact pressures reach their peak centrally and diminish toward the periphery.^17,19 In Carter and Wong's finite element model, the maximum load was represented as a distribution of pressures acting over a localized region of the cartilage surface, based on their consistency with those typically observed in adult joint contact mechanics. Although they did not report the total applied force explicitly, we estimated it by averaging the pressures applied to the elements. The loaded region spans a total length of 1.62 mm across 12 elements, yielding a loading radius (ri) of 0.81 mm. The sum of there reported pressures on elements corresponds to 54.8 MPa, resulting in a mean pressure of 4.57 MPa per element. Given that the model is two-dimensional, the pressures are applied over a line, which requires assigning a representative thickness to extrapolate the total force in three dimensions. We selected a nominal thickness (T) of 1 mm, consistent with conventions in similar modeling studies, which permits converting the line-integrated pressure into an equivalent total force (Fmax) of 7.4 N (Equation (12)). Importantly, the selected thickness does not influence the computed distributed load, since it is used both to convert the pressure field into a total force and subsequently canceled out during the calculation of ki (Equation (11)), which involves dividing by thickness. Therefore, thickness serves only to scale the magnitude of the applied load without affecting the applied mechanical stimulus.(12)Fmax=2riT∑j=112pj12=2⋅0.81mm⋅1mm⋅54.8MPa12=7.4N

To describe the magnitude of the remaining pressures, we introduce a new parameter: the amplitude of the maximum load kmax. The intensities of all other loads are defined as percentages of this peak value, as illustrated in Fig. 10. This parameter establishes a link between the load magnitude and its parabolic representation. A higher kmax leads to a more pronounced parabola, with a larger enclosed area, corresponding to a higher load. In essence, kmax serves as a quantitative link between the intensity of the mechanical stimulus and its parabolic representation. By applying the expression derived for ki (Equation (12)), an expression for this new parameter is obtained.(13)kmax=3Fi4riT.

Close

Fig. 10 Schematic representation of the loading history distribution for the Single and Double Contact Cases. Each parabola corresponds to an individual load defined by parameters ri, hi, and ki, scaled relative to kmax.

Export to PPT

Additionally, the loading history is defined as a sequence of such loads, each associated with the parameters ri, hi, and ki introduced earlier (see Fig. 9). The first two parameters (ri and hi) are fully determined by the geometry: the location and extent of the loads depend directly on the shape of the articular surface. For example, for Single Contact Case, if the surface is concave, the loads are concentrated within the concavity; if convex, they are spread laterally (Fig. 8). The loading histories for the Single and Double Contact Cases are defined as follows. To address this, the magnitudes of the applied pressures are specified as fractions of the maximum value, and they are symmetrical about the vertical axis. In the Single Contact Case, the loads are evenly distributed along the contact surface, and the span of the loads is determined by a new parameter ℓ that represents the maximum extent of the loading history. For the Double Contact Case, there are two clusters of loads, with their centers separated by another variable, denoted as λ. All the parabolic load profiles have the same radius, ri (Fig. 10).

These variables, which govern the extent of the loading history, are dependent on the shape of the head, specifically its concavity, quantified by the parameter αh. Expressions for they were obtained by fixing a value for a concave case (αh ≤ −15°), a planar case (αh = 0) or a convex case (αh ≥ 15°) and interpolating over the intermediate values, yielding to the following expressions.(14)ℓ(αh)=Lmin+Lmid−Lmin5,if αh≤−15°,Lmin,if αh=0,Lmid+35(Lmax−Lmid),if αh≥15°(15)λ(αh)=Lmin+2Lmid3,if αh≤−15°2Lmid+Lmax5,if αh≥15°

The radius of the applied pressure distribution is defined separately for single and double contact cases, as follows:(16)rl(αh)=λ(αh)2,if αh≤−15°λ(αh)4,if αh≥15°(17)rl(αh)=Lmax−Lmid3.

Each individual load is applied separately and later, the OI for each point the domain is calculated applying its expression (Equation. (5)). Ultimately, the geometry and boundary conditions are defined.

2.6

2.6 Convergence analysis

Ensuring the accuracy of the FEM solution requires careful mesh setup. We conduct a convergence analysis by varying mesh density to evaluate the solution's stability. This helps identify a suitable mesh that balances accuracy and computational efficiency. In our case, we selected two geometries and tested five mesh sizes. As mesh density increases, we analyze the behavior of the ossification index (OI) distribution within the cartilage, determining the mesh density that yields accurate and reliable results (see Fig. 11). The model exhibits a singularity near the surface interface on the side, where the values do not converge. This region is excluded from the analysis, as the overall distribution of the OI in the rest of the computational domain approaches to a solution. However, a balance between computational resource usage and solution accuracy is needed. Then, from case ne = 16 to ne = 32, the overall result exhibit initial stabilization. Therefore, selecting a value between these two cases is suitable. For this reason, ne = 20 was chosen, resulting in 202 edges and 7360 faces. Appendix A.3.1 presents the graphical overview of the applied meshes and the corresponding convergence results.

Close

Fig. 11 Convergence analysis results. Appendix A.3.1 provides a complete visualization of the convergence analysis results.

Export to PPT

3 Results

Once the mesh convergence was validated, we analyzed the geometries obtained and the distribution of the ossification index (OI) across different geometries and loading configurations. High OI values in cartilage regions suggest a higher likelihood of secondary ossification center (SOC) formation. We show 5 representative results out of the 25 simulations run for each scenario, as they capture the main trends observed. The full set of 25 cases for each scenario can be found in the Appendix, section A.3.2.

Scenario 1, which includes heads ranging from convex to concave with varying cartilage thicknesses (Fig. 12a), shows as the head transitions to a higher concavity, the OI maximum shifts closer to the surface, particularly for higher cartilage thicknesses (tc = 3.5 mm), where a greater shift in the OI position is observed. Both increased convexity and greater cartilage thickness are associated with higher OI values. The results also suggest that ossification at the articular surface increases in more concave heads, as OI values on the surface continue to rise with growing concavity. This may be explained by the progressive shift of the SOC toward the surface, which reduces the influence of compressive stresses associated with joint contact.

Close

Fig. 12 Ossification patterns under varying geometric and loading conditions corresponding to single contact case. Rows show the geometry case (Convex to concave, concave, convex head) and columns show variations of geometric variables of these cases. See Table 1 to identify the values of the variables. If a parameter has only one value on the table, it remains constant across all columns. They also correspond to the diagonal values of the results matrices on the appendix. Refer to A.3.2 for a complete visualization of the analysis results.

Export to PPT

In the double contact case (Fig. 12b), thin cartilage and convex heads do not generate a clear SOC pattern. As tc increases, OI peaks begin to emerge near the ossification front. However, increased concavity inhibits this pattern, causing ossification to occur more peripherally near the condyles. A high OI region between the condyles is observed for αh = −15°. Overall, higher tc and convexity promote SOC formation, while planar and concave heads (αh < 5°) result in SOCs appearing at the base of the specimen. Additionally, in this case, articular cartilage appears across all configurations, as indicated by consistently low OI values along the contact surfaces.

We now turn to scenario 2, where the horizontal and vertical radii are varied in fixed concave heads, namely, where we have αh = −15° fixed, and rx and ry are varied. In the single contact case (Fig. 12a), SOCs consistently appear. Wider heads (larger rx) shift the SOC closer to the surface, and lower ry enhances OI intensity. In contrast, under double contact (Fig. 12b), two local maxima near the condyles appear in almost every case, except for low values of rx and high values of ry. Additionally, epiphyses exhibiting sharper condyles present higher ossification on the articular surface.

Moving to Scenario 3, similar to the previous case, but with a convex head (αh = 15°). Under single contact (Fig. 12a), a centralized and slightly elongated SOC appears (visually observed as an elliptical high-OI region centered below the contact zone), especially in narrower heads. Increasing rx reduces overall OI, flattening the ossification pattern. With double contact (Fig. 12b), the model predicts a downward shift in SOC location and an increase in OI near the ossification front, producing a disk-like pattern, based on visual inspection of the OI distribution maps. These results show that convex heads are more sensitive to changes in the width of the head than in its height, which suggests that the horizontal shape has a greater effect on SOC formation than its height.

4 Discussion

The present model enables the creation, simulation, and comparison of various endochondral ossification cases. It provides significant insights into bone ossification and predicts the formation of ossification centers, whether single or double, within the same framework. We applied NURBS curves to gain more flexibility in defining the geometry, allowing us to explore a wide range of shapes. To gain insight into the ossification pattern, we used a simplified approach using a 2D linear elastic model. Additionally, we designed two types of contact configurations to reflect the variety of synovial joints found in the body.

Results emulates previous research, where diarthrodial joints were investigated in general, and the results were found to be particularly consistent with the metacarpophalangeal joint anatomy, in terms of the location and shape of the SOC. Additionally, the low OI values on the surface and at the ossification front suggest the formation of articular and growth cartilage, aligning with those previous findings,¹⁷ where concave and convex morphologies produced centralized ossification zones. Conversely, other models⁵ emphasize the role of biochemical factors, which are not captured in the current framework.

Moreover, our model effectively emulates synovial joints. For example, in the knee, a concave head with two contact loads resembles the distal femur load scenario, which is subjected to loads from the meniscus and the patella (Fig. 13a). Our model predicts a single ossification center positioned away from the contact surface and closer to the ossification front (Fig. 12b, Scenario 1, column b), showing a single ossification front in this area and resembling clinical data (Fig. 13b). In contrast, the contact surface of the tibia is nearly horizontal, and our model reflects this characteristic with an intermediate case where the head angle αh is close to 0°, resulting in a single contact scenario (Fig. 12a, Scenario 1, column c).

Close

Fig. 13 Visual representations of knee joint anatomy: (a) Diagram showing the structure of the knee joint (original image from ⁵⁰), (b) MRI scan highlighting the internal features of the knee, showing developed ossification centers on the tibia, femur and patella (original image from ⁵¹).

Export to PPT

In the case of the elbow, which is a compound joint, in its anteroposterior plane view the humerus shows a shape comparable to the concave head used in our model, with loading conditions similar to the double contact configuration (Fig. 12b, Scenario 2, column d). Although the model assumes geometric and loading symmetry, which contrasts with anatomical asymmetry in joints like the elbow, yet the model still predicts the formation of multiple SOCs, such as those associated with the capitellum and the trochlea, corresponding to contact points with the radius and ulna (see Fig. 14a).

Close

Fig. 14 Visual representations of elbow joint anatomy: (a) Diagram showing ossification centers on elbow joint (original image from ⁵²), (b) Radiography highlighting the important structures on the elbow (original image from ⁵³).

Export to PPT

Moreover, a minority of geometries and loads are not able to represent anatomical structures or may yield results that do not coincide with biological observations, which are excluded from the analysis. For example, non-biological patterns, such as ossification on the epiphyses’ surface or large ossification near model singularities, may emerge. For example, when a concave head has a high rx and low ry, the resulting geometry can present a high ossification index on the top surface, not aligning with the general prediction of SOCs being embedded in the cartilage. Nevertheless, simplified mechanobiological models can still provide valuable insights into how bone formation occurs, while higher complexity models can be applied to explore skeletal pathologies.¹² Moreover, as pointed out in previous research, fluid flow is significant at the boundaries of cartilage, and considering fluid exudation on cartilage surfaces could significantly impact areas with high fluid flow.³¹ Other improvements could include the incorporation of molecular reaction-diffusion models, allowing for better analysis of biochemical impacts,^{5,14,20,21,41} and dynamic simulations to emulate growth.^5,18

In conclusion, we have developed a computational framework able to analyze a wide range of bone anlage shapes and loading conditions under a static approach. Unlike earlier models that explored parametric development with limitations in geometric flexibility, this model allows for the investigation of an even broader set of geometries and loading cases. In future research, it can be applied to dynamic models where shape evolution is considered. This could involve altering initial conditions or defining the evolution model in terms of the specified NURBS parameters or including additional ones. While the current study employs a static approach, future research could enhance the realism of the model by incorporating extra considerations, such as biochemical factors, cellular effects, fluid flow dynamics, dynamic loading conditions, and ossification.

CRediT authorship contribution statement

Cristian Rodrigo Bustamante-Porras: Conceptualization, Methodology, Formal analysis, Investigation, Validation, Visualization, Writing – original draft. Kalenia Marquez-Florez: Supervision, Writing – review & editing. Carlos Alberto Duque-Daza: Supervision, Writing – review & editing. Diego Alexander Garzón-Alvarado: Conceptualization, Methodology, Supervision, Project administration, Writing – review & editing.

Guardian/patient consent

Not applicable.

Ethical statement

Not applicable.

Funding statement

The authors received no specific funding for this work.