A Review on Recent Advances in the Constitutive Modeling of Bone Tissue

Purpose of Review Image-based finite element analysis (FEA) to predict and understand the biomechanical response has become an essential methodology in musculoskeletal research. An important part of such simulation models is the constitutive material model of which recent advances are summarized in this review. Recent Findings The review shows that existing models from other fields were introduced, such as cohesion zone (cortical bone) or phase-field models (trabecular bone). Some progress has been made in describing cortical bone involving physical mechanisms such as microcracks. Problems with validations at different length scales remain a problem. Summary The improvement of recent constitutive models is partially obscured by uncertainties that affect overall predictions, such as image quality and calibration or boundary conditions. Nevertheless, in vivo CT-based FEA simulations based on a sophisticated constitutive behavior are a very valuable tool for clinical-related osteoporosis research.


Introduction
The Biological Material "Bone" Bone is a living material, organized in a fascinating microstructural hierarchy [12,37]. Besides its purpose as a mechanical support for the musculoskeletal system, it serves as a calcium reservoir [7], and a container for bone marrow. Bones in the human body consists of two distinct structures-the compact cortical shell and the spongy trabecular (or cancellous) region on the inside (Fig. 1). Cortical thickness ranges from several tenths of a millimeter to some millimeters in humans, while a trabecula is typically around 50-300 μm thick [37]. Individual trabeculae assemble to form foam-like structures with a typical bone volume fraction of 5-30% depending on anatomical location and age [48]. In cortical bone, Haversian and Volkmann channels, lacunae, and canaliculae lead to a cortical porosity of 4-17% depending on age [2,3,6].
The underlying microstructure of bone consists of lamellae of 3-6 μm thickness that arrange as osteons in compact bone and as lamellar packets in trabeculae [14,35,39]. Those lamellae are formed from bone fibrils-a highly aligned and organized staggered arrangement of bone's basic constituents: collagen type I molecules, hydroxyapatite mineral platelets, water, and a small fraction of non-collagenous proteins [10,50].
This complex hierarchical structure, alongside with the chemical composition, leads to the remarkable mechanical properties of bone that can be observed at different length scale levels. Pathological changes in nanostructural characteristics get manifest macroscopically, e.g., in osteogenesis imperfecta, where a collagen defect leads to brittle bones. This multi-scale interconnection of bone structures makes the constitutive modeling challenging.

Mechanical Behavior of Bone
The response of biological tissues to mechanical loads has been the focus of biomechanical investigations for more than a century. On the material level, this response is described as the relationship between the size-independent measures of mechanical stress and strain. For a certain strain stimulus that This article is part of the Topical Collection on Biomechanics is imposed, the material is reacting with stress response and vice versa. This response is intrinsic to the material and generally depends on the magnitude, the rate, the direction, the duration, and the number of repetitions of the imposed stimulus.
The description of an apparent behavior of an inhomogeneous material like bone requires an a priori definition of the considered length scale (nano, micro, macro). Microstructural features, like individual fibrils or lamellae, behave differently than bone at the tissue or organ level. For the sake of conciseness, this article will focus on the bone at higher length scales. One could imagine a representative cubical volume element (RVE) with a few millimeters in size, containing either cortical bone including vascular pores or a subsection of the trabecular network denoted as trabecular bone (Fig. 1b, d). Bone tissue properties at this macroscopic length scale are the most relevant in musculoskeletal research like implant technologies, tumors, trauma, and osteoporosis. Even at this scale, the description of the trabecular network relies on the underlying properties of the trabecular tissue, which is itself considered in an RVE as shown in Fig. 1c.
When loading macroscopic bone specimens in a testing frame, they react elastically if the load does not exceed the yield limit. While there are no differences between tensile and compressive elastic stiffness for compression and tension [38], there are differences in yield stress, strength, and fracture strain, Fig. 2a [5,13,34]. Due to its microstructure, bone's mechanical properties are direction dependent (i.e., anisotropic) but can be approximated with transverse isotropic or orthotropic material symmetry [9,20,33]. Consequently, when testing specimens from long bones in a longitudinal direction, a higher stiffness, yield stress, and strength can be observed, compared with testing in the transverse direction (Fig. 2a). Determining all 5 (transverse isotropy) or 9 (orthotropic) elasticity constants of bone experimentally is a very difficult task, for which ultrasonic measurements or a combination of tensile/torsion/shear tests can be used [54,55].
When loading bone past the yield point, two energy dissipation mechanisms occur: plasticity and damage Fig. 2b. Plasticity is manifest as permanent deformation after unloading, whereas damage is associated with a reduction of stiffness due to the formation of microcracks [21,44,47]. In contrast, the appearance of macroscopic cracks is denoted as a fracture. At very high compressive strains, trabeculae collapse and densification appear which leads to rising stresses.
Bone mechanical behavior is influenced by viscosity and is therefore time dependent. More precisely, its stiffness, yield limit, and strength depend on the strain rate. These quantities are larger, the higher the speed at which load is applied. At the same time, at very high strain rates, the bone becomes less ductile and the fracture strain is reduced (Fig. 2c) [8,18,25]. Moreover, bone also exhibits creep-the accumulation of plastic deformation under constant load-provided that the load exceeds a certain threshold (Fig. 2d).
These mechanisms of elasticity, plasticity, damage, viscosity, and creep are constitutive mechanical properties of the material bone and intrinsic in nature. They are always present but appear in different accentuations depending on the hydration state, the molecular composition, and the microstructure of the respective bone specimen. For example, dry bone tissue is more brittle and less viscous than bone with physiological water content. Furthermore, it is well known that the apparent stiffness, yield stress, and strength are reduced when porosity increases. This is evident for pathologic bone alterations like osteoporosis or cancer, in which the bone mass is reduced and subsequently, porosity is increased compared with healthy bone [4,15,48,52]. (the fraction of actual bone volume in a RVE) of trabecular bone specimens and stiffness as well as yield stress. Both properties are severely reduced for osteoporotic and cancerous BV/TV values [27]. Although not clearly visible in Fig. 3, usually, a nonlinear relationship exists between material parameters and BV/TV. It is known as a "power law density-elasticity relationship" that has usually the form  [24] and used with permission from Elsevier. b Uniaxial compressive cyclic behavior of trabecular bone showing yielding, stiffness reduction due to damage, and the remaining plastic deformation after unloading. The figure was adapted from [21] and used with permission from John Wiley and Sons. c The dependency of stiffness, yield, and failure on strain rate was obtained from compression tests on human cortical bone. The figure was adapted from [25] and used with permission from the American Physiological Society. d Creep behavior of cortical human bone. Three stress regimes are shown. For low stress, no creep occurs. For stress above a certain threshold, plastic creep strain is accumulating over time. The figure was adapted from [11,41] and used with permission from Elsevier a b Fig. 3 Dependency of stiffness a and yield stress b of human trabecular bone on BV/TV (bone volume/tissue volume). Healthy, osteoporotic, and cancerous BV/TV regimes are indicated. Adapted from [27], with permission from Springer Nature with E being the elastic modulus and A and B being constants [56][57][58].

The Principle of Constitutive Modeling
A constitutive material model is aiming to replicate what has been observed in reality by means of a mathematical framework based on fundamental physical principles (Fig. 4). In addition to replicating, it should be predictive, which means that although being derived from a cloud of discrete sample data, the model should be able to calculate responses that have not been tested experimentally. This allows for its use for, e.g., implant design or clinical evaluation of bone strength.
In the sense of mimicking a material response to a mechanical stimulus, a constitutive material model is operating in the scheme of stress and strain. It attempts to return a stress output similar to the real bone material when subjecting it to a certain strain input stimulus (or vice versa) (Fig. 4).
A constitutive model as such generally captures certain classes of constitutive behavior, e.g., elasticity, plasticity, and damage as described in "The Biological Material "Bone"" section.
To get a good match with the individual material, a number of model parameters must be adapted accordingly (Fig. 4b). Those parameters can be (i) material parameters like Young's modulus or yield stress, (ii) microstructural parameters like BV/TV or degree of anisotropy, or (iii) parameters with no physical meaning like exponents in power laws.
Before applying a new constitutive model to a practical problem, it needs to be properly validated. By validation, it is confirmed that the model is performing as expected and its output is in acceptable accordance with real-world observations. Usually, model validation is performed either in a dedicated set of experiments or with respect to another already validated model. Due to experimental restrictions, a validation might not be carried out on the length scale of the model but on a higher length scale, which is introducing additional uncertainties.
Considering the complex structure and behavior of bone tissue from "The Biological Material "Bone"" and "Mechanical Behavior of Bone" sections, constitutive material models for bone can be classified according to (a) length scale (RVE size, i.e., considered structure like cortical bone, trabecular bone, single trabeculae, homogeneous, or heterogeneous), (b) material symmetry (e.g., isotropic, orthotropic), (c) loading capabilities (small or large strain amplitudes, monotonic or cyclic amplitudes, low or high strain rates), and (d) constitutive response (elastic, plastic, viscous, damage, fracture, densification).
Constitutive models that cover a wide predictive range of material behavior are difficult to formulate mathematically and even more difficult to validate experimentally. Thus, constitutive models are formulated for a certain range of applications in which they describe reality within an accepted accuracy. If used outside those limits, the model will return data which is then deviating from what would be observed in an experiment more than can be accepted. For example, a simple linear elastic material model should not be used for strains that exceed the elastic limit of the material.

Review of Recent Constitutive Models
The second part of this review deals with the question if improved or new constitutive models have been developed and to what extent these findings influence bone strength prediction or related problems in musculoskeletal research.

Literature Search
A systematic literature search was done on PubMed. The intention was to review recent advances with respect to constitutive bone models. Thus, the search was limited to the past 5 years. An overview of constitutive models prior to 2016 can a b Fig. 4 a Principle of constitutive modeling. Material behavior is observed in reality by subjecting the material to an input stimulus (e.g., strain) and observing the output (e.g., stress). By creating a constitutive model, it is attempted to obtain a similar output for the same input and-at the same time-obtain trustworthy outputs for new input stimuli within a certain range of applicability. Parameters are used to tune the model to a certain type of material. b Stress-strain curves illustrating an exemplary experimentally observed behavior σ observed and the calculated model behavior σ calculated . Here, an elasto-plastic model is shown with a fair match of the observed data be found in [40•]. Studies were excluded if they did not present a reasonable new or improved constitutive model for bone.

Results
Selected publications on constitutive modeling of various types of bone are summarized in Table 1 and associated validations (experimental or numerical) are presented in Table 2.
In the following, individual aspects will be reviewed in more detail.
Predictive Range of Models Most models include plasticity, damage, and/or fracture. Viscosity and densification are rarely implemented. Usually, the models cover a specific range of predictable constitutive behavior. Elastic Limit/Plasticity A common challenge is to define the onset of elasticity in the case of multi-axial stresses. Yield surfaces based on Drucker-Prager, von Mises, or principle strain are typically used. Specifically, in the reviewed works, the elasticity limit is mostly based on a quadratic criterion of Drucker-Prager type and the post-yield behavior includes linear isotropic hardening. Some authors applied theories from other engineering areas such as the Lee-Fenves model from civil engineering [28], phase-field models from material engineering [45], or cohesive zone models from composite engineering [26]. Only the latter study provided a reasonable high number of experiments and showed fracture load correlations of R 2 = 0.89 for the presented cohesive zone model. In general, the multi-axial experimental validation of yield surfaces requires bi-or triaxial testing equipment and is thus very challenging [F, G].
Micro-architecture On the apparent trabecular bone level (Fig. 1b), architectural features like BV/TV and anisotropy (DA) are still the gold standard independent parameters to predict apparent stiffness and strength based on these CTbased measures [31,32] and two study approaches were  Fig. 1b; tb-ti, trabecular tissue Fig. 1c; ct-bo, cortical bone Fig. 1d; iso, isotropic; ortho, orthotropic; trans, transverse iso; N/A, not applicable; DPY(Q), Drucker-Prager yield, quadric approx. [42]; DPY, Drucker-Prager yield; DLY, Drucker-Lode yield; PSD, principle strain damage [29]; LFD, Lee-Fenves plastic damage [22]; 1DY, one-dimensional yield; BM, Burger model; VMY, von Mises yield; VM+C+D, von Mises yield+cap+ element deletion; CIY, cast iron yield; BD, Brittle damage [17]; DPY(Q)+D, DPY(Q)+densification [19]; CZM, cohesive zone model; PFM, phase-field model; DPD+D, DP+Q damage onset+element deletion; GMM, generalized Maxwell model found: first, direct homogenization of the apparent level by using bone density (power laws [16,45]) and, in some cases, the orientation of the trabecular network (fabric [31,32]). Power laws are simple to implement and are the choice if only bone density is available, whereas orthotropic fabric-based models are more accurate but such models need information about the trabecular bone orientation. In the second case, trabecular bone micro-FE models were used together with local constitutive modeling at the tissue level. However, such modeling requires high-resolution imaging as input and a direct validation was not done. Cortical bone was considered as homogeneous, isotropic, or orthotropic material and microstructures like pores are not considered explicitly. Damage and Fracture Interestingly, models including fracture mechanics became increasingly popular in the last years. Either rather simple "element deletion" models [46,51] or more sophisticated fracture energy-based models [26,45] have been published. Such approaches are more and more used in engineering fields and are available in commercial software packages which makes their use relatively easy. However, local "crack" modeling strategies include a process zone or mesh dependency at high strain amplitudes. This introduces a further-sometimes invisible-tuning parameter. Put simply, the in-elastic modeling effort is shifted from the development of a sophisticated constitutive model (including such processes) to a numerically more expensive simulation model which gives a good insight into the progression of fracture processes.
Anatomical Samples Material laws based on experimental studies provide parameters that depend on the biological material used but also on the shape of the sample and other factors like the number of samples and applied loading. Many studies used human femora bone but also studies on animal bones (bovine, porcine) were published for experimental comparison. Some studies were based on existing experimental results. As usual in the case of constitutive developments, many experiments (Table 2) were carried out on trabecular and cortical cubes/cylinders with a size of approx. 3-20 mm or whole bone parts. Only a few studies presented a sufficiently high number of specimens to capture the anatomical variations.

Model Validation
To begin with, it should be mentioned that outcomes presented in Table 2 which are based on a few specimens should be rather seen as a proof-of-concept than a validation. A direct multi-axial validation of clinically applicable apparent bone models is currently not possible experimentally. Therefore, as mentioned above, CT-based micro-FE models are used together with computational homogenization to predict the apparent behavior. However, Baumann et al. [1] stated that it is unlikely that a tissue constitutive model can be fully validated from apparent level experiments alone, as calculations are too insensitive to identify differences in the outcomes. In this context, the results of Sabet et al. [40•] demonstrated that using three different tissue constitutive models had only a slight effect on the apparent response. As expected, there was a significant change in the apparent response which goes back to the bone volume fraction. In the case of viscosity, Ojanen et al. [30] found that bone tissue visco-elasticity may not fully be able to predict the macro-scale visco-elasticity in human trabecular bone.
Software In order to use the developed constitutive laws for real bone models, they were implemented in finite element codes. Although many software packages enable nonlinear material modeling, only a few different packages are reported. Mainly, the commercial software Abaqus (Dassault Systemes Simulia Corp, Providence, RI) is used. It allows for any nonlinearity and the possibility of user material modeling. Parallelizable codes like Feap (University of California, Berkeley) and ParOSol (ETH Zürich, Switzerland) also formed the basis for some studies which provide even more implementation flexibility and a better parallel execution performance [46].
Clinical Applicability Most of the presented constitutive models can be used for bone strength assessment based on quantitative information about bone density. However, applicable new or sophisticated engineering approaches with respect to constitutive modeling of bone were not identified.

Discussion
The mechanical response of bone depends on factors like the structure or length scale under consideration and the applied loading regime. Constitutive models try to describe this response through mathematical models. Usually, such models are implemented into finite element tools to allow the analysis of complex geometries, material distributions, and loading conditions. The focus here was to evaluate recently published constitutive laws for trabecular or cortical bone with respect to practical usage, e.g., in CT-based bone strength predictions.
In summary, there are advances in the field of constitutive modeling of bone tissue. The main focus of recent studies was on the extension and experimental validation of existing models. There were a few new developments, but no fundamental improvements were found compared with the years before. Most important seems the consideration of physical mechanisms (e.g., microcracking) in the case of combined loading. In this regard, a generalization to 3D and application on bone structures which are subjected to cyclic overloading, e.g., bone implant systems, should be researched. Interestingly, a trend towards fracture-based models is observed over the last years. Such models put less effort into the local constitutive modeling but could lead to more insight into the overall failure mechanisms of bone. This might inspire new sophisticated constitutive model developments.
Some studies showed that a multi-scale validation (tissue to apparent level) is still inaccurate. Obviously, important features are missing or experimental shortcomings influence the results. For this reason, validations should take place on the same length scale on which the model is operating. For investigations at any level, more experimental data are desirable taking into account local strain measurements, multi-axial loading, and cyclic loading at different strain rates.
Model validations were carried out in vitro and showed similar accuracies as previous studies. One difficulty in recognizing the actual improvement through a new constitutive model lies in the fact that poor image quality and calibration or incorrect or oversimplified loading have a strong impact on the quality of the obtained results. This means there will be plenty of scope for further improvement of constitutive bone models, but at the same time, other uncertain model parameters must be reduced.
A general problem with simulation methods in biomechanics is the lack of standardization. Sample preparations, experimental protocols, tissue types, scanning, image processing, FEA modeling, and post-processing are done differently from study to study. It is not possible to compare constitutive models quantitatively making it impossible to separate good from less good models.
In terms of clinical applicability in osteoporosis research, in vivo CT-based FEA simulations based on sophisticated constitutive behavior still have great potential. Although no great achievements regarding constitutive modeling have been achieved, the FEA technology itself has developed accordingly. There is an increased level of automation and computational capabilities and the insights and predictive power of FEA simulations models are far better than estimates based on bone mass [49] or other non-physical image-based measures.
Funding Open access funding provided by TU Wien (TUW).

Compliance with Ethical Standards
Conflict of Interest Dieter H. Pahr is CEO and owner of Dr. Pahr Ingenieurs e.U. Pfaffsẗatten, Austria which develops and distributes the software medtool (www.medtool.at). Andreas G. Reisinger declares that he has no conflict of interest.
Human and Animal Rights and Informed Consent This article does not contain any studies with human or animal subjects performed by any of the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.