Simulating Metaphyseal Fracture Healing in the Distal Radius

Simulating diaphyseal fracture healing via numerical models has been investigated for a long time. It is apparent from in vivo studies that metaphyseal fracture healing should follow similar biomechanical rules although the speed and healing pattern might differ. To investigate this hypothesis, a pre-existing, well-established diaphyseal fracture healing model was extended to study metaphyseal bone healing. Clinical data of distal radius fractures were compared to corresponding geometrically patient-specific fracture healing simulations. The numerical model, was able to predict a realistic fracture healing process in a wide variety of radius geometries. Endochondral and mainly intramembranous ossification was predicted in the fractured area without callus formation. The model, therefore, appears appropriate to study metaphyseal bone healing under differing mechanical conditions and metaphyseal fractures in different bones and fracture types. Nevertheless, the outlined model was conducted in a simplified rotational symmetric case. Further studies may extend the model to a three-dimensional representation to investigate complex fracture shapes. This will help to optimize clinical treatments of radial fractures, medical implant design and foster biomechanical research in metaphyseal fracture healing.


Introduction
Accounting for 17% of all emergency department visits [1], distal radius fractures (DRFs) are the most common long-bone fracture. Its incidence is increasing because of an active lifestyle and an ageing society [1][2][3][4].
Typically, a DRF is a high-impact trauma that occurs in young patients during a fall, while in older patients it is characterized by a low-energy fracture occurring above the distal articular surface of the radius.
Particularly the tapering of the cortical bone, where the bone is reinforced by the trabecular structure, was found to be a critical point of failure [2]. DRFs can be categorized according to different classification schemes, depending on the type of defect and trauma [2,5,6]. DRFs can be treated either conservatively or by surgical intervention where the main objectives are repositioning, fracture reduction and immobilization to maintain the reduction [7,8].
The healing process of trabecular bone fractures differs from that of most diaphyseal fractures, in which healing occurs through endochondral ossification within a callus enclosing the fracture. In diaphyseal fractures, soft tissue is replaced by cartilage, which subsequently ossifies. In trabecular bone fractures, bone forms directly through intramembranous ossification with almost no intermediary cartilage [9,10]. However, in vivo studies by Claes et al. [11,12] and Sandberg et al. [13,14] demonstrated that cancellous bone healing follows similar biomechanical rules as diaphyseal healing even though the observed healing patterns differ.
The development of numerical models able to simulate such complex biological and mechanical processes has progressed significantly over recent decades [15][16][17][18][19][20][21][22][23][24]. But to date, most models have generally focused solely on the modelling of diaphyseal fractures in long bones. Shefelbine et al. [25] and Simon et al. [26] initiated one of the only attempts at modelling trabecular healing, later presented as the Ulm tissue-level bone healing model [21]. Thus, they predicted the sequence of events in metaphyseal fracture healing in an academic example with the use of a model able to predict diaphyseal and metaphyseal fracture healing.
As a part of a multinational collaboration over the course of the last few years, one of the most extensive clinical patient dataset tracking of in vivo fracture healing of the distal radius was acquired. Patients were inspected at six time points within the first year after fracture [27]. This allowed us for the first time to compare in vivo fracture data with the predictions of our fracture healing simulation.
Therefore, the aim of this study is to develop and adapt a pre-existing tissue-level bone-healing model to the metaphyseal fracture healing of the distal radius. To determine whether the Claes et al. [12] hypothesis that metaphyseal bone healing displays similar mechanobiological effects as diaphyseal healing appears reasonable because the developed simulation tool with its origin in simulating diaphyseal fractures provides an ideal possibility. Once we have established a reliable model, it may be useful as a tool for different research questions or even to provide support for optimizing clinical procedures.

Tissue Differentiation
The previously validated Ulm tissue-level bone healing model [21,25,[28][29][30][31][32][33] underlies the biomechanical observations of Pauwels [34] and a further mechano-regulating hypothesis by Claes and Heigele [35]. From mechanical quantifications, Simon et al. [21] derived local distortional and dilatational strains as the primary mechanical stimuli for the formation and removal of the following tissue types: woven and lamellar bone, cartilage and connective tissue. In addition to the mechanics, biological stimuli involved in the bone healing process are represented. Properties like the locally enclosing tissue types and level of vascularization augment the differentiation rules, which are evaluated by a fuzzy logic controller. Through the iterative numerical procedure, as depicted in Figure 1, the tissue differentiation processes of fracture healing in an arbitrary setting can be predicted. More detailed insights can be found in the Supplementary Materials to this paper and publications by Simon et al. [21] and Niemeyer et al. [28].

Figure 1.
The Ulm tissue-level bone healing model. The numerical procedure of solving an initial value problem to describe the biological tissue concentrations in an arbitrary setting uses a fuzzy logic controller. This process uses the input of a structurally finite element analysis, the enclosing tissue types and further biological settings to determine the concentration changes of woven and lamellar bone, cartilage and connective tissue. Further, the extensions of this study, the fractured trabecular (crushed trabecular) bone is depicted. The numerical procedure is adapted from Niemeyer et al. [18] CC-BY. The tissue differentiation diagram is adapted with permission of Georg Thieme Verlag KG (Stuttgart).

Metaphyseal Healing
As Claes et al. [12] demonstrated, both metaphyseal and diaphyseal bone healing appear to be regulated by the same biomechanical stimuli (dilatational and distortional strain). They found that metaphyseal areas, exposed to low strains led to intramembranous bone formation, whereas higher strains additionally provoked endochondral ossification or fibrocartilage formation. As this coincides with assumptions of the diaphysis, all current theories and tissue differentiation processes as implemented in the most recent healing model [28] were used and not adapted or further calibrated.

Homogenized Material for Trabecular Structure
The present model assumes a continuum, that is, organ-level, representation of the involved materials. Therefore, by modelling a metaphyseal bone healing process, the complex trabecular geometrical structure is not required, but its mechanical influence can be approximated by the bone volume content of trabecular bone tissue over time. Studies by Arias-Moreno et al. [36], Hosseini et al. [37], MacNeil and Boyd [38] and Varga et al. [39] evaluated the differences between a detailed, high-resolution finite-element model of the trabecular structure (µFE) in the metaphyseal radius and a representative homogenized material finite-element model (hFE). In those studies, strain analysis revealed that a homogenization of the trabecular structure is a sufficient approximation. Considering that representing all of the fine trabeculae as in a µFE analysis would increase the computational effort for a fracture-healing simulation to an insoluble problem. By this approach, the trabecular structure is represented by a spatially homogenized mixture of lamellar bone and soft tissue. The Ulm tissue-level bone healing model. The numerical procedure of solving an initial value problem to describe the biological tissue concentrations in an arbitrary setting uses a fuzzy logic controller. This process uses the input of a structurally finite element analysis, the enclosing tissue types and further biological settings to determine the concentration changes of woven and lamellar bone, cartilage and connective tissue. Further, the extensions of this study, the fractured trabecular (crushed trabecular) bone is depicted. The numerical procedure is adapted from Niemeyer et al. [18] CC-BY. The tissue differentiation diagram is adapted with permission of Georg Thieme Verlag KG (Stuttgart).

Metaphyseal Healing
As Claes et al. [12] demonstrated, both metaphyseal and diaphyseal bone healing appear to be regulated by the same biomechanical stimuli (dilatational and distortional strain). They found that metaphyseal areas, exposed to low strains led to intramembranous bone formation, whereas higher strains additionally provoked endochondral ossification or fibrocartilage formation. As this coincides with assumptions of the diaphysis, all current theories and tissue differentiation processes as implemented in the most recent healing model [28] were used and not adapted or further calibrated.

Homogenized Material for Trabecular Structure
The present model assumes a continuum, that is, organ-level, representation of the involved materials. Therefore, by modelling a metaphyseal bone healing process, the complex trabecular geometrical structure is not required, but its mechanical influence can be approximated by the bone volume content of trabecular bone tissue over time. Studies by Arias-Moreno et al. [36], Hosseini et al. [37], MacNeil and Boyd [38] and Varga et al. [39] evaluated the differences between a detailed, high-resolution finite-element model of the trabecular structure (µFE) in the metaphyseal radius and a representative homogenized material finite-element model (hFE). In those studies, strain analysis revealed that a homogenization of the trabecular structure is a sufficient approximation. Considering that representing all of the fine trabeculae as in a µFE analysis would increase the computational effort for a fracture-healing simulation to an insoluble problem. By this approach, the trabecular structure is represented by a spatially homogenized mixture of lamellar bone and soft tissue.

Rule of Mixture
Because the elastic modulus scales according to the cube of the mineral density [40], we developed a rule of mixtures based on the experimental findings of Carter and Hayes [40], as a result of which the averaged material properties of the tissue compositions in each point could be calculated. As a linear calculation of the Young's modulus is not appropriate, as demonstrated by Simon et al. [21], a cubic rule of mixture should be used. Therefore, the tissue composition of N-sorted tissues with elasticity parameters E 1 < E 2 < . . . < E N and the corresponding concentrations c 1 ...c N are combined following Pietsch et al. [41]:

Modelling of a Compression Fracture
Regarding the modelling of a compression fracture in a metaphyseal bone, the fractured area must be defined accurately. Therefore, we introduced here a new tissue type that we term "crushed trabecular bone" (CTB). As the name suggests, it should represent the fractured trabeculae bone structures, which are damaged and pressed into each other. In vivo and in vitro observations [11,14,[42][43][44][45][46][47] showed that this damaged structure still contributes to the primary stability of the bone but is resorbed over time.
Modelling this densification is achieved by increasing the bone concentration in the fractured area according to an exponential function (2). This approach follows in vivo observations and studies by Thurner et al. [48] and Hambli [49], who investigated the compression of trabecular structures. The calculated concentration was assigned to the new designed material type CTB: Thus, the fractured area consists of a certain amount of CTB and connective tissue. Material properties with a Young's modulus of 6000 MPa and a Poisson's ratio of 0.3 were used for pure CTB [50]. The densification approach, as well as the material properties, followed in vivo observations. Beside resorption, this new tissue type does not participate in any other tissue differentiation processes as depicted in Figure 1.
Niemeyer et al. [28] implemented a method to consider the influence of adjacent tissue on the local differentiation process. This distance-weighted average tissue potential (µ) was used to formulate an influence of the surrounding substances on the resorption speed of CTB (r CTB ). From this, each tissue type was represented by a specific resorption potential p ∈ [0, 1]. Tissues with high perfusion tended to contribute more to the resorption of the fractured bone structure. Consequently, higher resorption potentials (r) were assigned according to Table 1.
Here, the Gaussian kernel G adj is defined as introduced by Niemeyer et al. [28]. Resorption rates followed recent findings on the resorption speed by Atkins et al. [51] and Christen et al. [52]. Through parameter variations (see Section 3.1 Parameter Variation) using a different patient's geometry than the ones reported in Figure 4 and using comparisons with in vivo observations, it was concluded that a maximal resorption speed of 11% per day (r ctb = 0.11) was used in all simulations.

Application to Distal Radius Fractures
Because we have access to datasets of distal radius fractured patients, the newly developed method was applied to four patient specific geometries. Patients were scanned at six time points with a Scanco medical Xtreme CT-II-scan system in weeks 1, 3, 5, 12, 24 and 52 after trauma. All patients were treated conservatively and the cast was removed before the third scan. During the healing period, patients were urged to perform a daily hand gripping exercise using a hand-grip fitness tool. Approval was obtained from the ethical committee of the Medical University Innsbruck for this study (AN2014-0374 344/4.31).
To simulate this DRF treatment, a patient-specific simulation geometry, initial tissue concentrations and boundary conditions were defined according to Figure 3B.
Regarding computational efforts, two-dimensional (2D) rotational symmetric assumptions were used. Comparing the simplified 2D rotational symmetric simulation domain in the following patient-specific geometry application, it was important to note that the entire radius is represented by the 2D geometry. Only the radial styloid process cannot be modelled by this 2D approach ( Figure 2). Inaccuracies in the articular surface are covered by the principle of Saint Venant. The geometrical patient specification was achieved through parameterized geometry generation. Parameters such as cortical thickness, radial diameters and fracture site were measured in the scans after training and assistance by clinicians. Details on the investigated patients are attached in the Supplementary Materials.

Application to Distal Radius Fractures
Because we have access to datasets of distal radius fractured patients, the newly developed method was applied to four patient specific geometries. Patients were scanned at six time points with a Scanco medical Xtreme CT-II-scan system in weeks 1, 3, 5, 12, 24 and 52 after trauma. All patients were treated conservatively and the cast was removed before the third scan. During the healing period, patients were urged to perform a daily hand gripping exercise using a hand-grip fitness tool. Approval was obtained from the ethical committee of the Medical University Innsbruck for this study (AN2014-0374 344/4.31).
To simulate this DRF treatment, a patient-specific simulation geometry, initial tissue concentrations and boundary conditions were defined according to Figure 3 Regarding computational efforts, two-dimensional (2D) rotational symmetric assumptions were used. Comparing the simplified 2D rotational symmetric simulation domain in the following patient-specific geometry application, it was important to note that the entire radius is represented by the 2D geometry. Only the radial styloid process cannot be modelled by this 2D approach (Figure 2). Inaccuracies in the articular surface are covered by the principle of Saint Venant. The geometrical patient specification was achieved through parameterized geometry generation. Parameters such as cortical thickness, radial diameters and fracture site were measured in the scans after training and assistance by clinicians. Details on the investigated patients are attached in the supplementary material.  The Geometries are discretised by a finite element mesh with quadratic hexahedral elements (3724-4143 elements), and the mesh was checked by mesh convergence analysis.
Boundary conditions were chosen to approximate the real situation within the constraints of a 2D rotational symmetric simulation: The bone was fixed to the ground and compressed onto the articular surface ( Figure 3A). This artificial load case was intended to represent the daily mechanical stimulus to which the radius is exposed.
The Geometries are discretised by a finite element mesh with quadratic hexahedral elements (3724-4143 elements), and the mesh was checked by mesh convergence analysis. Boundary conditions were chosen to approximate the real situation within the constraints of a 2D rotational symmetric simulation: The bone was fixed to the ground and compressed onto the articular surface (Figure 3-A). This artificial load case was intended to represent the daily mechanical stimulus to which the radius is exposed. During the healing period, the patients used a hand-grip fitness tool (approx. 100 N compression force). The resulting axial compression force from the contracting muscles in the wrist served as the mechanical stimulus for the fracture healing process. The occurring force onto the radius was calculated with the aid of a recently developed musculoskeletal inverse dynamics model of the human hand [53][54][55]. Results showed that, on average, the radial joint surface was exposed to an axial load of 351 N without great influence of the patients' anatomical variability.
The initial model had homogenous material properties: trabecular bone, cortical bone and marrow. To obtain a realistic baseline tissue distribution within the trabecular bone, the model underwent tissue differentiation until a homeostatic situation was obtained (Figure 3-C). In this case, the apparent stiffness of the tissue concentrations was sufficient, and the implemented differentiation processes achieved no change in tissue concentrations. Through this process, the homeostatic situation of the specifically determined geometry and boundary condition was obtained as exemplary depicted in Figure 3-C.
Through X-ray scans, an approximation of the patients fractured area was identified and assigned a CTB concentration as defined by formula (2) and depicted in Figure 3-D in blue.
To simulate the conservative treatment until week 5, a rotational symmetric brace with specific material parameters [56] (Young's modulus of 950 MPa and Poisson's ratio of 0.3) was modelled around the radius geometry (Position (a) in Figure 4). Muscles, skin and other soft tissue in between are assumed to be pure connective tissue. To compare the simulation with six in vivo longitudinal X-ray scans, the simulation results were visualized by an emulated pseudo-X-ray view. Concentration values of lamellar and woven bone and CTB were multiplied by factors for their corresponding bone mineral density [57].
The entire simulation framework was implemented in Python (Version 3.7, Python Software Foundation, Beaverton, OR, USA) using the commercial finite-element solver ANSYS (Release 2020 R1; ANSYS, Inc., Canonsburg, PA, USA). Running the entire simulation procedure for one geometry on a local, well equipped workstation results in calculation times of approximately 180 min. During the healing period, the patients used a hand-grip fitness tool (approx. 100 N compression force). The resulting axial compression force from the contracting muscles in the wrist served as the mechanical stimulus for the fracture healing process. The occurring force onto the radius was calculated with the aid of a recently developed musculoskeletal inverse dynamics model of the human hand [53][54][55]. Results showed that, on average, the radial joint surface was exposed to an axial load of 351 N without great influence of the patients' anatomical variability.
The initial model had homogenous material properties: trabecular bone, cortical bone and marrow. To obtain a realistic baseline tissue distribution within the trabecular bone, the model underwent tissue differentiation until a homeostatic situation was obtained ( Figure 3C). In this case, the apparent stiffness of the tissue concentrations was sufficient, and the implemented differentiation processes achieved no change in tissue concentrations. Through this process, the homeostatic situation of the specifically determined geometry and boundary condition was obtained as exemplary depicted in Figure 3C.
Through X-ray scans, an approximation of the patients fractured area was identified and assigned a CTB concentration as defined by formula (2) and depicted in Figure 3D in blue.
To simulate the conservative treatment until week 5, a rotational symmetric brace with specific material parameters [56] (Young's modulus of 950 MPa and Poisson's ratio of 0.3) was modelled around the radius geometry (Position (a) in Figure 4). Muscles, skin and other soft tissue in between are assumed to be pure connective tissue. To compare the simulation with six in vivo longitudinal X-ray scans, the simulation results were visualized by an emulated pseudo-X-ray view. Concentration values of lamellar and woven bone and CTB were multiplied by factors for their corresponding bone mineral density [57].  The entire simulation framework was implemented in Python (Version 3.7, Python Software Foundation, Beaverton, OR, USA) using the commercial finite-element solver AN-SYS (Release 2020 R1; ANSYS, Inc., Canonsburg, PA, USA). Running the entire simulation procedure for one geometry on a local, well equipped workstation results in calculation times of approximately 180 min.

Results
The predicted fracture healing results of four patient-specific distal radius fractures are compared in Figure 4 by an emulated pseudo-X-ray view to longitudinal in vivo X-ray scans of the corresponding patient.
Clearly, artefacts of the conservatively treated patients' braces were visible in the first two screenings before the cast was removed prior to the third scan. In the fractured areas, the trabecular structure was compressed. During weeks 5 to 12, this compressed structure was resorbed and new bone formed. In most cases, the fractured area was totally healed after one year. In Patient A, the higher bone density in the area of the compression fracture clearly decreased from day 21 to day 81. Simultaneously, new bone was formed until day 168. The simulated results displayed comparable results. The compression fracture, modelled through an increased concentration of CTB, was resorbed in the first 3 to 5 weeks and was simultaneously replaced by new bone.
In Position (a), the modelled brace around the radius was depicted and was present until week 5. Patient A at day 36 (Position b) displayed newly formed woven bone in the fractured area, but this bone was not evenly distributed. As the CTB structure was resorbed and laterally a higher concentration of CTB was still present, not as much woven bone was formed in this spot, but even here the woven bone formation concentration attained 53% compared to a maximal 86.4% in the medial area. This observation was also obtained in the corresponding X-ray, where specific spots of the fractured area were also more mineralized than the surrounding. When comparing the subjects to each other, it was noted that larger fractured areas tended to heal more slowly. In addition, the predicted homeostatic tissue distribution showed that smaller radii tended to exhibit higher bone density than larger radii.
The time course of tissue concentrations in the fractured area revealed the underlying healing processes in more detail: Figure 5 demonstrates the mean concentration values in the fractured compartment of Patient D. First, the CTB was rapidly resorbed over time. After 59 days, the resorption speed decreased because from days 29 to 65 woven bone was formed. This mostly originated from intramembranous ossification, but because only little cartilage formation was observed at the beginning some endochondral ossification also occurred. Clearly visible was the maturation process, whereby the newly formed woven bone became lamellar bone, which was, after 178 days, the only tissue presents apart from connective tissue.

Parameter variation
Formula (3) introduced a new parameter into the differentiation processes: the maximal resorption speed of CTB ( ). Variations of this parameter of between 5 and 23 % (representing tissue loss per day) are depicted in Figure 6. A higher resorption rate resulted in drastically decreased bone volume. Bone formation initially increased in all cases and subsequently decayed to the same concentration level.

Parameter Variation
Formula (3) introduced a new parameter into the differentiation processes: the maximal resorption speed of CTB (r ctb ). Variations of this parameter of between 5 and 23% (representing tissue loss per day) are depicted in Figure 6. A higher resorption rate resulted in drastically decreased bone volume. Bone formation initially increased in all cases and subsequently decayed to the same concentration level.

Parameter variation
Formula (3) introduced a new parameter into the differentiation processes: the maximal resorption speed of CTB ( ). Variations of this parameter of between 5 and 23 % (representing tissue loss per day) are depicted in Figure 6. A higher resorption rate resulted in drastically decreased bone volume. Bone formation initially increased in all cases and subsequently decayed to the same concentration level.

Discussion
In this study, a novel model to simulate metaphyseal bone healing based on the tissue differentiation hypothesis of Claes and Heigele [35] was introduced. The underlying Ulm fracture healing model for diaphyseal fracture healing [21,28] was used and extended for its application to the metaphysis. In addition to geometrical and boundary adaptations, the model had to be substantially extended to consider further materials and processes to simulate the observed in vivo effects.
The newly introduced material, CTB, represents the compressed and fractured trabecular tissue. Trabeculae which are broken and of which the blood supply is interrupted, cannot contribute to the accumulation of mesenchymal stem cells and the resulting bone formation processes. CTB comprises a passive component. Because its material properties are similar to those of trabecular bone, it stabilizes the fractured area and allows rapid bone formation by its mere existence. Furthermore, the modelled CTB had no direct biological influence, as proposed by Aspenberg and Sandberg [13], in that no obvious relationship between newly formed bone and the old trabeculae was observed. Despite the fact that biopsies [13] showed that osteoconduction on old trabeculae is not an important process, we maintained the consideration that its surface offers ideal osteoconductive properties [58].
The modelled CTB was resorbed through a single rate in all biological processes. This rate was dependent on the surrounding tissue type, which might be an appropriate assumption, but modelling more specific resorption processes could be promising.
Modelling the geometry of the radius as a rotationally symmetric representation on the one hand might be a good approach for limiting computational costs, but on the other hand it limits the ability to represent a patient-specific radius shape. The radial styloid process as well as the asymmetrical shape of the lunate and scaphoid facet cannot be considered in such an assumption. On top, the fractured area and its positioning are limited to a horizontal expansion and can only be personalized through its height and vertical positioning. However, using 3D computations would have increased the computational time by up to 2 weeks for each healing simulation, but that is unacceptable in the development process of a new method.
Nevertheless, the model also shows that in 2D the assumptions on metaphyseal fracture healing appear reasonable, and further studies with 3D models can be elaborated using the same framework and modelling techniques.
From a large dataset, the four presented patients ( Figure 4) were selected because they covered a broad variety of different aspects. The geometrical shape as well as their sex and age varied widely. The selected patients' X-ray scans were compared to the simulation outcomes through an emulated X-ray presentation. Therefore, it must be noted that the emulated image was closer to a computed-tomography scanned image slice than to an actual X-ray view because the 3D geometric radioscopy was not modelled in this view. We nevertheless believe that a comparison of the real and simulated results is clearer through such a presentation. Because of the long time-lag between X-ray scans, a detailed comparison of the healing process was difficult. A qualitative comparison showed the model to be capable of simulating the metaphyseal section, the compression fracture and its healing progress. The simulations were able to predict the healing patterns observed in the X-ray scans from a decrease in bone density in the fractured area to new bone formation over the same time frame.
Discrepancies in the bone concentration of Patients A and D around the rotational axis (Positions c and d) might have arisen from the artificial load case as well as the boundary conditions required to maintain axial rotation symmetry. The significant difference in the diameter between Patients B and D was reflected in the bone concentration (Positions e and f) in the proximal region and might be explained by the fact that the model first calculated a homeostatic situation based on hand gripping force, which triggered bone formation in the relatively small cortical shell of Patient C (see (e) in Figure 4).
De Jong et al. [43] observed an increase of total bone (11%) and trabecular bone (20%) density until weeks 6-8. In the simulated results, a similar behaviour was observed. Increases of total bone concentration of 11 to 17% from days 10 to 90 were predicted in our simulations.
Published quantified results of in vivo investigations on distal radius fractures by de Jong et al. [42,43] were compared to congeneric quantities obtained by the simulations (Figure 7). Same behaviours were obtained for the axial stiffness, with a decrease from day 21 to 42. Density comparisons showed that the increase of healthy woven and later lamellar bone slightly lagged behind the in vivo observations but showed the same behaviour.
De Jong et al. [43] observed an increase of total bone (11 %) and trabecular bone (20 %) density until weeks 6-8. In the simulated results, a similar behaviour was observed. Increases of total bone concentration of 11 to 17 % from days 10 to 90 were predicted in our simulations.
Published quantified results of in vivo investigations on distal radius fractures by de Jong et al. [42,43] were compared to congeneric quantities obtained by the simulations (Figure 7). Same behaviours were obtained for the axial stiffness, with a decrease from day 21 to 42. Density comparisons showed that the increase of healthy woven and later lamellar bone slightly lagged behind the in vivo observations but showed the same behaviour. Histological studies by Aspenberg and Sandberg [13], Chen et al. [47], Jarry and Uhthoff [9] and Claes et al. [12] showed that metaphyseal bone formation occurs mainly by two patterns. Either lamellar bone forms on existing trabeculae, or woven bone emerges directly in the marrow [13] between the old trabeculae without direct connection to preexisting bone [47]. Cartilage was seldom observed [12,13], but islands of chondrogenic tissue indicated limited endochondral ossification [12]. All observations showed that direct intramembranous bone formation is the dominant process of bone formation.
Those findings strongly correspond to our simulated results. We detected little cartilage formation in the fractured area ( Figure 5) and hence little endochondral ossification. Direct woven-bone formation in the fractured area indicated that the mechanical stimuli were within the ranges originally proposed for the diaphysis. The course of tissue formation in Figure 4 thus reflected the histological observations that new woven bone formed in between the old trabeculae (CTB) and matured within a plausible time to lamellar bone. Additionally, the amount of bone formed corresponded to the surrounding bone concentration of the uninvolved parts of the bone. X-ray scans one year after trauma showed a healed radius, and it was difficult to determine the previously fractured area. The same characteristics occurred in the simulated results.
The histological comparison showed that the simulation model represents the metaphyseal fracture healing in an accurate manner. Because the model assumes tissue differentiation processes and biomechanical findings originally formulated for diaphyseal fracture healing, it can be stated that this simulation model corroborated the hypothesis of Claes et al. [12]. Despite the fact that some parameters of the model might differ slightly between the diaphysis and metaphysis, the same mechanobiological rules controlling bone and cartilage formation and resorption are able to predict fracture healing patterns both for the diaphysis and the metaphysis. Histological studies by Aspenberg and Sandberg [13], Chen et al. [47], Jarry and Uhthoff [9] and Claes et al. [12] showed that metaphyseal bone formation occurs mainly by two patterns. Either lamellar bone forms on existing trabeculae, or woven bone emerges directly in the marrow [13] between the old trabeculae without direct connection to preexisting bone [47]. Cartilage was seldom observed [12,13], but islands of chondrogenic tissue indicated limited endochondral ossification [12]. All observations showed that direct intramembranous bone formation is the dominant process of bone formation.
Those findings strongly correspond to our simulated results. We detected little cartilage formation in the fractured area ( Figure 5) and hence little endochondral ossification. Direct woven-bone formation in the fractured area indicated that the mechanical stimuli were within the ranges originally proposed for the diaphysis. The course of tissue formation in Figure 4 thus reflected the histological observations that new woven bone formed in between the old trabeculae (CTB) and matured within a plausible time to lamellar bone. Additionally, the amount of bone formed corresponded to the surrounding bone concentration of the uninvolved parts of the bone. X-ray scans one year after trauma showed a healed radius, and it was difficult to determine the previously fractured area. The same characteristics occurred in the simulated results.
The histological comparison showed that the simulation model represents the metaphyseal fracture healing in an accurate manner. Because the model assumes tissue differentiation processes and biomechanical findings originally formulated for diaphyseal fracture healing, it can be stated that this simulation model corroborated the hypothesis of Claes et al. [12]. Despite the fact that some parameters of the model might differ slightly between the diaphysis and metaphysis, the same mechanobiological rules controlling bone and cartilage formation and resorption are able to predict fracture healing patterns both for the diaphysis and the metaphysis.

Conclusions
This metaphyseal bone-healing model provided defined and adjustable biomechanical conditions for a wide range of radius geometry representations. The simulations confirm the hypothesis of Claes et al. [12] for similar mechanobiological rules in fracture healing of both, the metaphysis and diaphysis. Furthermore, the results demonstrated reasonable predictions of intramembranous, and little endochondral, ossification in the fractured area. Therefore, the model appears appropriate for studying fracture healing under differing mechanical conditions and varying metaphyseal bones and fracture types. Investigating comorbidities such as osteoporosis and the implicated reduced bone density, as well as the extension of the model to a 3D representation where patient-specific, nonsymmetric fractured areas can be investigated, are part of future studies. Those fracture healing simulations could then help to understand the healing pattern better and derive clinical treatments as well as biomechanical knowledge of metaphyseal fracture healing. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.