A comparative study of orthotropic and isotropic bone adaptation in the femur

Functional adaptation of the femur has been studied extensively by embedding remodelling algorithms in finite element models, with bone commonly assumed to have isotropic material properties for computational efficiency. However, isotropy is insufficient in predicting the directionality of bone's observed microstructure. A novel iterative orthotropic 3D adaptation algorithm is proposed and applied to a finite element model of the whole femur. Bone was modelled as an optimised strain-driven adaptive continuum with local orthotropic symmetry. Each element's material orientations were aligned with the local principal stress directions and their corresponding directional Young's moduli updated proportionally to the associated strain stimuli. The converged predicted density distributions for a coronal section of the whole femur were qualitatively and quantitatively compared with the results obtained by the commonly used isotropic approach to bone adaptation and with ex vivo imaging data. The orthotropic assumption was shown to improve the prediction of bone density distribution when compared with the more commonly used isotropic approach, whilst producing lower comparative mass, structurally optimised models. It was also shown that the orthotropic approach can provide additional directional information on the material properties distributions for the whole femur, an advantage over isotropic bone adaptation. Orthotropic bone models can help in improving research areas in biomechanics where local structure and mechanical properties are of key importance, such as fracture prediction and implant assessment. © 2014 The Authors. International Journal for Numerical Methods in Biomedical Engineering published by John Wiley & Sons, Ltd.


INTRODUCTION
Adaptation algorithms have been incorporated into finite element (FE) studies in many areas of biomechanics that focus on bone morphogenesis and response to altered loading conditions [1][2][3][4][5][6][7]. Bone was initially assumed to be a self-optimising linearly elastic continuum that responded to changes in strain energy density (SED) [1,6,8,9]. Coelho et al. [10] and Kowalczyk [11] have used SED as the driving stimulus for the optimisation of bone with a hierarchical macrostructural and microstructural description. However, SED can produce convergence problems during the adaptation process at a continuum level [4]. The action of directional-dependent normal strains on the bone matrix has been put forward as the generator of physiological mechanobiological signals that 874 D. M. GERALDES AND A. T. M. PHILLIPS activate osteocytes [12,13] and better suited as the driving stimuli of the adaptation process in continuum models [4,14,15]. In order to model the process of bone adaptation, the driving stimulus needs to be a physiologically meaningful representation of the in vivo mechanical environment [16]. Therefore, the FE model of the bone being studied is required to be as close to the physiological state as is reasonably possible. This involves careful selection of its constitutive representation, mesh, geometry, loading and boundary conditions [17]. FE models of the femur frequently fall short of these criteria. In order to simplify the analysis, 2D representations of the femur [11,15,18,19] and partial models [6,10,14,16,[20][21][22][23] are commonly used and ignore the adaptation process in different planes or regions of importance. Artificial boundary conditions, such as restricting displacement at a distal end of the femoral shaft, induce stress concentrations around the restrained region [5, 6, 14-16, 20, 22-24]. Non-physiological loading conditions such as applying hip contact forces (HCFs) and muscle forces as point loads [10,14,16,21,25,26] are often adopted for simplicity. These affect the strain and stress distributions in the surroundings of the load application point as well as in the trabecular bone beneath the cortex, so increased discretisation has been recommended [27]. Phillips [28] proposed a free boundary condition approach to produce an equilibrated model of the femur, with the ligaments and muscles spanning across the hip and knee joints explicitly included. Balanced models such as this do not directly constrain the bone in a non-physiological manner [27,29] and produce reduced absolute strain values for the femur undergoing a single-leg stance [24,27], resulting in more physiological stress and strain distributions [30].
Bone is usually modelled with isotropic material properties in an attempt to reduce computational times [1,2,31], despite the anisotropic nature of the material properties being measured experimentally [32][33][34]. Orthotropy has been shown to be the closest approximation to the bone's anisotropy, short of full anisotropic modelling [32]. In addition, isotropy is insufficient in predicting the directionality of the observed microstructure of the bone [35][36][37][38].
The need for a physiological continuum model of the material properties distribution and structure orientation across the femur in order to understand its biomechanical behaviour has been emphasised [39]. A review of the regression equations that have been fitted between elastic properties measured experimentally and computed tomography (CT) derived densities suggests that it is difficult to accurately determine this relationship [40]. Furthermore, CT images are composed of scalar density values resulting from a combination of local porosity and tissue mineralisation and, therefore, are not able to predict the directionally dependent elastic properties of the bone required to model its structural directionality at a continuum level [41]. Recent developments in micromechanics and X-ray physics [42] have allowed for extraction of orthotropic elastic properties from CT data. These studies rely on observer-dependent estimations of the trajectories of the principal material directions from the bone's geometry and from recognisable collagen structures amongst volumetric CT data of varying resolution [43][44][45].
Consequently, an iterative strain-driven orthotropic 3D bone adaptation algorithm was developed as part of the presented study. Orthotropic orientations of each element were aligned with the local principal stress directions, following Wolff 's Law [37,38,46]. Directional material properties were updated proportionally to the local strain stimuli, according to the Mechanostat theory for bone adaptation [47], which has been shown to have a logical biochemical basis [48][49][50].
The orthotropic adaptation algorithm proposed was compared against the commonly used isotropic approach. Both approaches were applied to a fully balanced loading configuration with the muscles and ligaments spread along their attachment sites. Boundary conditions that do not constrain the femur in non-physiological ways were adopted. The predicted density distributions for a coronal section of the whole femur were qualitatively and quantitatively compared with the results obtained by the commonly used isotropic approach to bone adaptation as well as ex vivo imaging data. The study investigates how important the contribution of orthotropic adaptation can be in the production of bone models where directionally dependent mechanical properties are required to be physiologically represented. We hypothesise the following: (1) the orthotropic approach will produce material distributions closer to that observed in vivo; and (2) additionally, orthotropy can produce reliable information about the local structural composition of the bone.

Finite element model
The FE model of the femur was built in Abaqus/CAE (v. 6.12, SIMULIA, Dassault Systèmes SA, Vèlizy-Villacoubly, France). The external geometry of the femur was extracted from CT scans of a third-generation Sawbones synthetic femur (Pacific Research Laboratories, Inc., Vashon, WA, USA) and meshed with 326 026 four-node linear tetrahedral elements. These had a mean element edge length of 2.37 mm. Mesh density and properties were within values found to be adequate for numerical modelling of the femur by a previous study [51]. Sensitivity studies were also performed by the authors to assess the dependence of the predicted results on mesh properties [52]. The femur was positioned with the centre of the femoral head coinciding with the origin of the global reference system and the x, y and´axes defined according to International Society of Biomechanics recommendations [53]. The femur's positioning and orientation were extracted from HIP98 [54], where kinetic and kinematic data were measured in vivo in combination with hip joint contact forces from instrumented prostheses. The peak frame of the fourth trial of the level walking load cycle for the subject HSR was selected as the load case to apply on the FE model (HSRNW4). The data produced from this subject's trials have been considered reliable by other studies [55,56]. In order to equilibrate the femur during the static analysis without influencing muscle recruitment, the inertial action of the contralateral leg and the missing torso was considered [55]. The inter-segmental moments and forces between the pelvis and the femur were extracted from the inverse dynamics analysis of the HIP98 kinematic data by an open-source musculoskeletal model validated against the same public dataset [55,57]. These were then applied at the centre of the hip joint structure and can be seen in Table I.
Some modifications were made to the femur model to promote the transfer of loads across the hip, greater trochanter, tibio-femoral and patello-femoral joints back to the femur. Bi-layered structures composed of a 1-mm-thick internal isotropic elastic layer representing cartilage (E D 5 MPa, D 0.49; Figure 1, in red) [58] and a 1-mm-thick external isotropic elastic cortical bone layer (E D 18 GPa, D 0.3; Figure 1, in grey) [7] were projected from hand-picked  contact surfaces for each of the three joints and modelled with linear six-node wedge elements. The external surface nodes of the bi-layered structures were connected to two nodes located medially and laterally along the functional axes for the tibio-femoral and patello-femoral joints ( Figure 1, in orange) or at the centre of rotation for the hip via sets of linear beam elements with elastic properties similar to the cortical bone and a cross-sectional area of 10 mm 2 . The functional axes used for the tibio-femoral and patello-femoral joints were measured by Klein Horsmann in a cadaveric study [59] and were introduced in order to restrict the joints' movements to the flexion/extension plane, in a hinge-like fashion, a common assumption made in musculoskeletal models [55]. At the condyles and the patella, the two nodes along the functional axis were connected via a beam element in order to transfer the moments between them ( Figure 1, in green). The medial condylar node was fixed against displacement, allowing the femur to pivot about the tibial plateau at the tibio-femoral joint [60] (Figure 1(a)). Stiff two-node axial connector elements were included to create a hinge-like trapezoid structure, representing the region of the patella connected to the patellar ligament and the quadriceps and allowing for force transfer back to the femur ( Figure 1(b)). At the hip joint, an artificial stiff truss structure connected the acetabular region to muscle insertion points on the pelvis, sacrum, lumbar spine and a point representative of L5S1, as proposed by Phillips [28] (Figure 1(c)). The hip was modelled as a pin joint with the inter-segmental forces and moments defined in Table I applied at its centre of rotation. A total of 26 muscles and seven ligamentous structures were explicitly included. Their origination areas as defined in the muscle standardised femur [61] were mapped onto the femoral mesh and the number of connectors, C , that formed each group was given by the number of nodes within these mapped areas. The insertion points in the pelvic and tibial regions were extracted from Phillips' [28] model of the femoral construct [28]. The typical force-displacement curve for these connectors was defined by a reference stiffness value, k ML iso , and a peak contractile force, F M peak , from [62] ( Figure 2). The values for k ML iso for each muscle group were calculated according to Equation (1), where L T slack is the tendon slack length [63].
Stiffness was lowered after 0.75F M peak , promoting activation of other muscle groups as the force generated by each muscle approaches its maximum [28]. Given that muscles act in tension [63], stiffness values were considered to be negligible (0.01k ML iso / in compression. The properties of the muscles and ligaments included in the model are compiled in Table II. Changes were applied to the definition of some muscles in order to represent their anatomical position and function more physiologically, in comparison with the model proposed by Phillips [28]. An extra connector element spanning between the turning point at the greater trochanter and the insertion point at the tibia was included in the iliotibial tract structure. Its stiffness was also increased to 97 N/mm, following tensile tests by Merican [64]. The peak force for the patella tendon was increased to the 2500 N measured as its ultimate failure load [65]. The tensor fascia latae was modified to allow for wrapping around the greater trochanter and force transfer back to the femur. This region of the muscle contact was modelled as a bi-layered structure, similar to the joint structures described previously (Figure 1(c)). Fixed constraints were applied at the insertion points of the muscles and ligaments on the tibia and fibula [28]. The point in space representing the insertion at the lumbar spine was connected to ground via a spring element with a stiffness of 10 N/mm in the anterior-posterior direction and negligible stiffness in the other two directions (0.1 N/mm), in order to simulate the balancing action provided by the upper torso during different gait activities [66]. The resulting model with all the joints, muscles and ligamentous structures can be seen in Figure 3.

Adaptation algorithms
The complete femur model presented earlier was loaded and submitted to one of two different optimisation processes: orthotropic or isotropic adaptation.

Orthotropic adaptation.
At each iteration, the stress, ij , and strain, " ij , tensors were extracted for each element and post-processed in MATLAB (v. R2007b, MathWorks, Natick, MA, USA). An eigenanalysis of the stress tensor gives the local principal stress vectors, pv , and corresponding principal stress values, p (Equation 2).
The element orthotropic material orientations were rotated to match with the calculated local principal stress orientations, following Wolff's trajectory theory [38,67] and previous work for 2D [15,52,68], in agreement with proven optimum orientations for orthotropic materials undergoing a 879 single load case [69]. The three local orthotropic axes, x 1 , x 2 and x 3 , were associated with the minimum, min pv ; medium, med pv ; and maximum, max pv , principal stress vectors, respectively. The local strain stimulus associated with the transformed orthotropic material axes, " ij , was calculated for each element according to Equation (3) [52,68].
The material properties were adjusted in order to bring the local strains within the remodelling plateau, as proposed in the Mechanostat theory for bone adaptation [47] A normal target strain, " nt , of 1250 strain with a margin of 0 D˙0.2" nt , was used to define this plateau. Each iteration, the Young's modulus, E it i , of the elements outside the remodelling plateau was updated proportionally to the absolute value of the associated normal local strain stimulus (Equation 4), and limited between 10 MPa and 30 GPa [41].
Shear moduli, G it ij , were taken to be a fraction of the mean orthotropic Young's moduli [15,52,68] Poisson's ratios for each element, it ij , were assumed to be less than or equal to 0.3 [70,71]. In order to satisfy the thermodynamic restrictions on the elastic constants of the bone, some elements' Poisson's ratios were altered, ensuring that the compliance matrix remained always positive definite (Equation 6) [67].
If E it j was greater than E it i , it ij was kept at 0.3, while it j i was adjusted such that the equality constraint was maintained (Equation 7).
A 'dead zone' of elements was defined in order to exclude elements of low elastic stiffness and undergoing negligible strains from the applied convergence criteria. An element was considered to be in the 'dead zone' when the maximum absolute normal strain values, " i i were less than 250 strain and, simultaneously, its directional Young's moduli, E n i , were all below 100 MPa. The adaptation process was considered to achieve a state of convergence when the average change in Young's moduli of all elements outside the 'dead zone' was less than 2% between successive iterations. All elements were assigned the same initial orthotropic elastic constants MPa/ and local orthotropic orientations matching the global femoral axis system. A convergence study was performed to confirm that the model had a limited range of sensitivity to this arbitrary starting point [52].

Isotropic adaptation.
Because isotropic symmetry does not require any information on direction of the material properties, Young's modulus for each element, E it , was updated proportionally to the maximum absolute value of the principal strains according to Equation (8). Similar to the orthotropic adaptation, the isotropic Young modulus values were updated to bring the local strains within the plateau of˙0.2" nt around the same normal target strain, " nt , and were limited between 10 MPa and 30 GPa. Based on the isotropic assumption, G it was taken as Equation (9).
Poisson's ratio for each element, it , was assumed to be equal to 0.3. A 'dead zone' was also defined, and the same convergence criteria for the orthotropic adaptation was applied. All elements were assigned the same initial isotropic constants .E D 3000 MPa, D 0.3, G D 1500 MPa/.

Imaging data
The results produced by both the orthotropic and isotropic adaptations were compared with ex vivo imaging data of two femur specimens. A CT scan of an ethically obtained specimen of a male cadaveric femur, 55 years old, weight 94 kg and height 188 cm, was taken with a SOMATOM Definition AS+ scanner (Siemens AG, Munich, Germany) based in the Queen Elizabeth Hospital, Birmingham, UK. The specimen was scanned at 120 kV and 38.0 mAs with an effective spatial resolution of 0.71 mm. The density and normalised density greyscales of the coronal slice of the whole femur were compared with the predictions for isotropic and orthotropic adaptations.
In addition, micro-CT (CT) data for the proximal region of a male cadaveric femur, 27 years old, weight 75 kg and height 175 cm, were also ethically obtained. The specimen was scanned using a HMXST 225 CT cone beam system with a 4MP PerkinElmer Detector (Nikon Metrology, Tring, UK) based in the Natural History Museum, London, UK, at 145 kV and 150 A and with an effective spatial resolution of 78.7 m. The directionality of the trabecular structures of a coronal section of the femoral head was compared with the dominant orthotropic orientations. Table III summarises the predicted values obtained for the resultant components (F r , F x , F y and F´/ of the HCFs for both converged models in %BW. The isotropic model produced higher resultant HCFs than the orthotropic model (332% against 319%). The isotropic model took three  fewer iterations to converge than the orthotropic model (26 vs 29). Both models achieved the predetermined convergence criteria, with an oscillating behaviour observed beyond 10 iterations as they attempted to reach the remodelling plateau.  For the converged solutions of the orthotropic and isotropic iterative simulations, the density of each element was calculated using a modulus-density relationship measured for the trabecular bone in the femoral neck [72]

RESULTS
The maximum relative density value achieved in the iteration process was 2.41 g/cm 3 . The average representative Young's modulus and density for both models are shown in Table IV. The converged isotropic model was 60% stiffer and 25.4% denser than the orthotropic model. The average values for the three main regions and the whole femur are also included.  Figure 4 shows the coronal slice representations of the density distributions resulting from the isotropic (a) and orthotropic (b) adaptation processes for the converged femur. All elements with density above 1.4 g/cm 3 were grouped together as dense cortical bone, in order to allow for better visualisation of the predicted density distributions for the trabecular bone. These were compared with a coronal slice taken from the CT scan (c) of the whole femur. Figure 5 shows the greyscale profiles of nine slices across the femoral head and neck region, shaft and condyles (transverse slices taken at 5%, 20%, 40%, 60%, 80% and 95% of the length of the femur) for the three coronal representations of the density distributions seen in Figure 4. These profiles were extracted using ImageJ [73], normalised and plotted between 0 and 1 (the maximum relative density value) along the percentage width of the slices taken.
The root mean squared error (RMSE, %) and Pearson's product-moment coefficient (r/ between the two different predictions and the CT scan profiles can be seen in Table V. These were calculated for the first third (medial aspect, 0-33%), the last third (lateral aspect, 66-100%) and the whole width (0-100%) of the slice. The average values for the three main regions and the whole femur are also processed. The converged orthotropic model predictions were in general found to be closer to the CT greyscale profile than the converged isotropic model. The orthotropic predictions resulted in lower RMSE and higher Pearson's r, particularly for the lateral side of the slice widths, indicating a better correlation with in vivo measurements. Poorer predictions (r < 0.5) were produced in the femoral condyle region for both models. Both material symmetry assumptions produced good results for the medial cortical shaft, whilst the lowest correlations between profiles were found in slices 1, 6 and 9, cutting through regions largely composed of trabecular bone.
For the orthotropic model, the dominant material directions were defined as the orientation associated with the highest directional Young modulus for each element. These were projected onto the density distributions for coronal representations of the proximal femur (Figure 6(a)) and compared with a 5-mm-thick reconstruction of CT slices of the same region (Figure 6(b)).

DISCUSSION
The loading environment of the FE model has been shown to be a key stepping stone in the attempt to produce a physiologically meaningful driving stimulus for the adaptation process [14,16]. Table III shows that both isotropic and orthotropic approaches produced relatively small differences between HCFs, agreeing with other comparisons between mechanical environments for FE models with isotropic and orthotropic material properties [74,75]. The HCFs predicted were higher than the values measured in vivo for the same instance of the gait cycle [54]. However, computationally derived HCFs are often over-predicted in comparison with measurements from instrumented prosthesis [28,55]. Figure 4 shows a comparison between the isotropic (a) and orthotropic (b) model predictions and a coronal slice of the CT scan of the whole femur specimen (c) for the same region. These predictions show good agreement with the structures produced by similar 2D [1,15,19] and 3D adaptation studies [20,23] of the proximal femur. Many important density features observed in the CT scans were correctly predicted by both approaches, such as the low-density regions in the intra-medullary canal, Ward's and Babcock's triangles, the region just above the principal compressive group in the femoral head and medial to the superior part of the femoral neck and the epicondylar regions for both condyles. However, the isotropic assumption did not accurately represent the density distribution for a coronal slice, as it did not predict the dense cortical distribution along the lateral shaft, superior aspect of the femoral neck and the greater trochanter. The orthotropic predictions of the thick femoral shaft and the denser regions in the superior and inferior regions of the femoral neck, along the surface of the greater trochanter and around the intercondylar notch are seen to be coherent with the CT scan. The use of a modulus-density relationship (Equation 10) to produce the density plots in Figure 4 is limiting because it is a specimen-dependent empirical relationship obtained for the femoral neck. Because this relationship is used for both the isotropic and orthotropic model results, the comparison between them still holds. More rigorous approaches have been put forward [43,44] and would allow for a more accurate representation of the density distribution in the femur.
The orthotropic assumption produced more defined Ward's and Babcock's triangles. It also resulted in a less-stiff and lower-mass optimised femur, under the same loading conditions (Table IV). Predictions of greyscale profiles for nine slices across the femoral head, shaft and condyles show that the orthotropic predictions have higher Pearson's r coefficients for the lateral shaft and the slices across the femoral neck. The RMSE value is generally lower for the profiles predicted with orthotropy compared with the isotropic assumption. However, in regions where mainly trabecular bone is present (such as slices 6 and 9, across the condyles and femoral head, respectively), the profiles show poor agreement compared with the CT scans. These are regions where trabeculae have been proposed to adapt to the shear stresses arising from complex loading scenarios [37,76]. The inclusion of a shear modulus adaptation algorithm could, therefore, have a positive impact in the density distribution predictions. It would also overcome the limitation of the ad hoc assumption of taking shear moduli to be a fraction of the mean orthotropic Young moduli.
A further advantage of using orthotropic material properties instead of isotropic symmetry is that directionality of the bone material properties can also be predicted. The proposed continuum approach we present in this study circumvents the assumption of using a pre-defined library of microstructure geometry [11] because it allows the system to optimise the combination of material orientations in order to provide the minimum energy solution for the load case it is subjected to. The adaptation algorithm proposed attempts to represent the underlying behaviour of both the cortical and trabecular bone structures as an orthotropic continuum with optimised material properties and orientations, whilst representing the complete femur with all muscle groups and ligaments explicitly modelled and load applicators included to promote physiological load application, with the possibility of being developed for multiple load cases in the future. Figure 6 shows the dominant orthotropic axes predicted by the algorithm (a) in comparison with a CT reconstruction for a similar coronal slice of the proximal and distal femurs (b). The directionality of all main trabecular groups documented was predicted for the proximal femur [36]: (i) the principal compressive group, composed of stiff, densely packed trabeculae arching from the medial cortex of the shaft towards the articular surface; (ii) the secondary compressive group arising from the medial cortex of the shaft, right below the principal compressive group, and curving towards the superior region of the femoral neck and the greater trochanter; (iii) the principal tensile group, composed by trabeculae arching from the lateral cortex across the neck of the femur and ending in the inferior part of the femoral 885 head; (iv) the secondary tensile group, originating in the lateral cortex just inferiorly to the principal tensile group, arching towards the mid-line of the upper end of the femur; and (v) the greater trochanter group, made by poorly defined trabeculae along the greater trochanter region. These arise as a structural response to the necessity to transfer load along the femur from an oblique to vertical direction [35,76,77]. Some other interesting features were also satisfactorily predicted. The cortical thickness seems to be a good match with the one observed in the CT slice, and trabeculae can be seen radiating from the centre of the femoral head and meeting its superior surface at perpendicular angles [35,78].
Despite the encouraging results, there are several limitations to the current models and methods. The geometry of the FE model of the femur was representative of neither the geometry of the subjects from which HCFs were measured nor the specimens used to obtain the CT and CT scans. Subject-specific hip geometry has been found to influence the calculation of the HCFs [79]. The definition of certain muscles as straight paths of action and the inter-segmental forces used could have contributed to the overprediction of the HCFs [55]. These reasons may explain some of the discrepancies between the calculated forces and the ones measured in vivo [54]. It can also give a justification for some of the differences observed in the predicted density distributions and trabecular directionality. The displacement restrictions implemented at the medial condyle node introduced artificial constraints to movement in the distal part of the equilibrated FE model, possibly resulting in estimation errors of the loading environment and, therefore, influencing the predictions for this region. The modelling of surface contact at the hip and knee joints could induce a more physiological behaviour of the model.
It is estimated that the average person undertakes 1.1 million walking cycles a year [80]. The load case selected corresponds to the peak of the level walking load cycle, where contact forces exceeding 250% body weight going through the femur have been measured in vivo [54]. This load case has been used in several FE and remodelling studies with good results for the proximal region of the femur [2,3,14,16,21,26] but may accentuate the differences between the isotropic and orthotropic models. The use of a single load case is a significant limitation when describing the complex mechanical environment the femur is subjected to physiologically [37,81], particularly in the distal region, also adapted to higher flexion load cases [76]. Reduced wall thickness in the anterior and posterior aspects of the cortical shaft may be a result of the simplified loading applied. Further work should include more load cases for a variety of frequent daily activities in order to improve the prediction of the distribution of the mechanical properties and associated orientations, particularly in the distal part of the femur and across the femoral head. These are regions that have evolved to resist the shear stresses arising from the multiple load cases the femur undergoes [35,37,76,78]. Inclusion of more daily activities may result in an increase in the shear stresses in the femur, with the resistance of the trabecular structures to shear potentially playing a more important role than in the single load case model [37]. Therefore, a shear modulus adaptation algorithm will be included in future studies. The proposed model results in an optimised structural system to most efficiently resist the loads generated by the instance of the walking cycle with highest contact forces, much like the trabecular structure that has been thoroughly studied in the proximal femur [35][36][37]. Time and frequency dependence of the orthotropic adaptation process would need to be considered if the model were to be extended to the study of bone remodelling around implants and bone fractures.
The topological density distributions predicted by the model seem to be in reasonable agreement with data extracted from the CT scan. The directionality generated at the continuum level was comparable with CT slices along the coronal plane for the proximal femur region. It also showed to be coherent with published studies of the trabecular features of these regions [35,36,38,78]. We believe the method proposed in this study can provide an alternative way of assigning orthotropic material properties to continuum models of bone, with material orientations aligned to resist the principal stresses arising from the loading the femur is subjected to.
It is tentatively concluded that the orthotropic assumption is more advantageous in comparison with the isotropic material symmetry assumption, confirming the initial hypothesis of the study outlined. Orthotropy provides a more accurate representation of bone's elastic symmetry and can also give information about the three-dimensional directionality of bone's tissue-level material properties. The use of a balanced model allows for the prediction of the adaptation process for the whole femur, without artefacts induced by the application of fixed boundary conditions directly on the bone in question. An orthotropic model for the complete 3D femur has been produced. The inclusion of multiple load cases and of a shear modulus adaptation algorithm could further improve the predictions. A robust orthotropic continuum model of the whole femur has potential in achieving a more thorough understanding of bone's structural material properties, thus improving the knowledge we have of its mechanical behaviour and response to the various loading environments it may be subjected to. Such a model could contribute to the improvement of the design of orthopaedic implants and fracture fixation devices, providing information on the directional properties of the bone surrounding these devices and how it may adapt to the changing mechanical environment.