Numerical homogenization of trabecular bone

The following report introduces various simulation campaigns aimed at homogenizing trabecular bone tissue.
Background, analyzing the structure of the trabecular bone it is evident how the shape of the pores and the trabecular structure influence the stiffness and the structural response. It is analyzed how the variation of the pore geometry and their distribution influences the homogenization parameters. In addition, specific case analyzes are carried out.
Results: a homogenization matrix and various analyzes on the variation of its parameters related to the variations of the pores are obtained. The most influential parameter is the fraction of pores but there can also be a large contribution given by the variation of the axial ratio;
Conclusions, the results obtained show how bone remodeling strongly contributes to increase the stiffness in the load direction. A homogenization matrix to consider both constituents together can be used but with the foresight to respect the volumetric fraction and verify that the axial ratio is close to 1.

The following report is extracted from a project of the Computational Mechanics of Tissues and Biomaterials course, held for Medical Engineering at the University of Rome Tor Vergata.

Introduction

The purpose of the following analysis is to investigate the influence of microstructural properties on the macroscopic behavior of the trabecular bone. The logical flow of the simulation campaigns is shown in fig. 2.

We start by analyzing the generic structure of the trabecular bone, arriving at choosing an idealized representative volume (RVE) with a geometry that respects the relationship between the pores and the trabecular bone mass. This RVE is then discretized and therefore the next step foresees to optimize the choice of the mesh in such a way as to have a sufficiently accurate result by reducing the computational effort requirements. A first homogenization matrix is ​​then obtained which reflects the structural properties of the RVE and trabecular bone thus idealized.

Different geometric and material properties are then analyzed following the responses of the bone to the load, a phenomenon known as bone adaptation and remodeling. In this simulation campaign it is analyzed how the variation of pores, their volume fraction and their distribution influences the macroscopic mechanical behavior.

Further analyzes are carried out to see how the axial relationship, the variation of the material parameters of the single elements, such as bone matrix and pore content, affect the structure. Further analyzes are then made by varying the window of interest and considering a specific patient case.

Bone tissue

Bone tissue is the biological tissue that forms bones. It is characterized by a remarkable hardness and resistance, due to the mineralization of the extracellular matrix. It is a constantly evolving fabric, continuously renewed and remodeled for the duration of its life. It constitutes the scaffolding of the skeleton, from insertion to muscles and tendons and contains bone marrow inside.

Bone tissue is made up of cells and an extracellular, organic and inorganic matrix. The organic part contains collagen fibers and amorphous substance formed by glycoproteins and proteoglycans.

Two types of bone can be distinguished: compact bone and cancellous bone [1]. The compact bone appears, on macroscopic examination, as a single solid mass. The spongy bone, on the other hand, has an alveolar appearance and is made up of thin trabeculae, a set of small plates and rods that anastomose into a three-dimensional network. The bone marrow is housed within the meshes that are formed with the trabecular spines. Trabecular bone constitutes 20% of the total bone mass and is located within the bone structures. Its distribution varies greatly from bone to bone.

Cortical and trabecular bone differ not only for their architecture, but also for the speed of development, function, proximity to the bone marrow, blood circulation, rapidity in turnover, fractures and response to stimuli. All of these factors provide different material behaviors. It should be emphasized that the bone, as a whole, consists of both a part of compact bone and a part of cancellous bone, so the response of the entire structure will have both contributions.

Spongy bone is typically less resistant than compact bone. Many tests have been conducted in the literature to evaluate the material properties of cancellous bone and there are many results, some even controversial. Typically, tests are performed on: compression, bending, uniaxial traction, ultrasonic measurements, backward calculation through 2D / 3D FEM simulations and microindentation [2]. The range for the elastic modulus is from 076 GPa to 20 GPa. It is now accepted that the trabecular bone module is 20-30% lower than that of compact bone.

Furthermore, it is known in the literature that the fundamental and most influential parameter in binding the microstructure to the stiffness response is the volumetric fraction of pores. [3]

But there are other things too. It must be taken into account that bone is a biological tissue and therefore subject to remodeling. In particular, it is known that bone responds very well to compression stimuli by increasing its cross section and thickening the compact bone layer. Some changes are also seen in the cancellous bone, including the alignment of the trabeculae along the lines of force. [4]

Multiscale structure

In this analysis, a lumbar vertebra, L3, is used as a reference. The considerations are however extensible to any other trabecular bone but the material parameters refer to data selected for lumbar vertebrae. The lumbar vertebrae are just some of the 33-34 vertebrae that make up the human spine. The vertebral body consists mainly of trabecular bone with the individual trabeculae aligned along the lines of force. This arrangement makes an important contribution to load distribution and torso support. Approximately three quarters of the entire structural load falls on the anterior column, that is, on the vertebral body, intervertebral discs and on the end plates of the vertebrae. The trabecular bone structure can undergo various changes with aging or due to certain pathologies [4] such as the loss of mineral density, morphological changes such as an increase in the volumetric fraction of pores, thickening or thinning of the trabeculae and loss of connectivity. These changes greatly affect the mechanical response of the tissue.

The bone is made up of different phases. It contains an organic phase, an inorganic phase and water. The organic phase is composed of collagen and proteins while the inorganic phase is calcium phosphate, similar to hydroxyapatite [5]. At the nanoscale, collagen molecules and hydroxyapatite crystals are found. We then move on to a slightly larger scale where mineralized fibrils unite to form a single lamella a few μm thick. The individual lamellae pack together to form trabeculae with an average thickness of 50μm and a length of approximately 1mm. We are on the microscale. We then move on to an intermediate scale where we can consider the volume of interest of this analysis where we find both trabeculae and pores filled with marrow up to the macroscale where we look at the bone as a whole, the individual material properties and the structure provide a well-defined structure response.
The goal of the following analysis is to take into account the differences at the microscale, considering a well-defined representative volume, and to homogenize the mechanical properties by obtaining a single homogenized stiffness matrix. This homogenized matrix describes the overall material behavior, taking into account both the individual material properties and their spatial and geometric distribution.

Homogenization

The goal is therefore to obtain material properties at the macroscale starting from the properties at the microscale such as the stiffness of the two main materials, bone matrix and pore content.

On a sufficiently large scale with respect to the individual trabeculae, it is possible to consider the pores, filled with bone marrow, and the bone trabeculae as two different materials but both heterogeneous and with isotropic behavior. The isotropy is a first approximation which however allows us to proceed more efficiently in homogenization by going to see how the geometry, the arrangement and the volumetric fraction of the pores affect the macroscopic behavior.

With the homogenization procedure we will consider an average of the mechanical properties of bone and matrix but this average is understood in a broader sense and the meaning may vary depending on the approach. The contributions of the individual material properties are influenced by the geometry of the structure and their distributions.

Volume of interest

The first step, in order to proceed with a homogenization, is to define a representative volume (RVE) such as to be representative from the statistical point of view of the geometric distribution of the constituents, therefore of the microscale. To do this, the code automatically calculates the side in the range of a few mm starting from some considerations. The mean pore diameter [6] and the volumetric fraction [2] are respected.

On the basis of these considerations, 25 equidistant pores are inserted and such that their distance from the edge is half the distance between each pore. This allows for a symmetrical and representative RVE of a selection of the trabeculae. More information on the procedure for automating the RVE calculation is in section 8.2.

In addition, to make the appropriate considerations in analyzing the average tensions, the pores are considered filled with a low-stiffness material representative of the bone marrow. This serves both to consider a phenomenon that is more representative of reality and to allow a correct application of the average stress theorem.

Analytical estimates

To compare the homogenization data, the analytical estimates of Voigt and Reuss are considered. The Voigt approximation considers the stiffness matrix as a weighted average for the volume fractions of the stiffnesses of the constituents.

\begin{equation}
\mathbb{C}_{\text {Voigt }}=\left(1-v_{f}\right) \mathbb{C}_{b}+v_{f} \mathbb{C}_{m}
\end{equation}

Reuss’s estimates are similar but compliant:

\begin{equation}
\mathrm{S}_{\text {Voigt }}=\left(1-v_{f}\right) \mathbb{S}_{b}+v_{f} \mathbb{S}_{m}
\end{equation}

These two estimates are considered due to Hill’s theorem which guarantees, under the appropriate hypotheses, that the stiffness matrix is included between them [7]. Reuss constitutes the lower limit and Voigt the upper limit of homogenization.

The two materials are considered individually isotropic, therefore such as to have a material behavior described by the compliance matrix in equation 3.

\begin{equation}
[\mathbb{S}]_{\text {iso }}=\left[\begin{array}{cccccc}
\frac{1}{E} & -\frac{v}{E} & -\frac{v}{E} & 0 & 0 & 0 \\
-\frac{v}{E} & \frac{1}{E} & -\frac{v}{E} & 0 & 0 & 0 \\
-\frac{v}{E} & -\frac{v}{E} & \frac{1}{E} & 0 & 0 & 0 \\
0 & 0 & 0 & \frac{2(1+v)}{E} & 0 & 0 \\
0 & 0 & 0 & 0 & \frac{2(1+v)}{E} & 0 \\
0 & 0 & 0 & 0 & 0 & \frac{2(1+v)}{E}
\end{array}\right]
\end{equation}

These matrices are calculated and for each iteration of the homogenization process it is verified that the result is included in these estimates using the VoigtReussTest [] function. In particular, passing the variable verbOutput == True it is possible to extract its values:

\begin{equation}
\mathbb{C}_{\text {Voigt }}=\left(\begin{array}{ccc}
2019.340 & 865.404 & 0 \\
865.404 & 2019.340 & 0 \\
0 & 0 & 576.969
\end{array}\right)
\end{equation}

\begin{equation}
\mathbb{C}_{\text {Reuss }}=\mathbb{S}_{\text {Reuss }}^{-1}=\left(\begin{array}{ccc}
0.2262 & 0.0399 & 0 \\
0.0399 & 0.2262 & 0 \\
0 & 0 & 0.09316
\end{array}\right)
\end{equation}

Numerical homogenization

A simulation campaign is then structured to analyze the RVE and extract a homogenized stiffness matrix. Further considerations are then made on the parameter C (om) 11. The homogenization approach follows homogeneous boundary conditions of the displacement type. In particular, a displacement field is applied at the macroscale assuming that it is representative of that at the microscale. Thus a certain voltage field is obtained at the microscale, a function of the point. Therefore it is possible to consider the voltage field at the macroscale as the average of that at the microscale. In particular, these conditions are guaranteed by the Average Strain Theorem. Therefore the mean stresses measured in the RVE are defined as:

\begin{equation}
\langle\sigma\rangle=\frac{1}{\|\tilde{\Omega}\|} \int_{\Omega} \varepsilon d \tilde{\Omega}
\end{equation}

where ˜Ω is the area of the RVE. By applying a deformation field to the macroscale, deformations are obtained at the microscale from which it is possible to extract the tensions simply by solving the problem of elastic equilibrium:

\begin{equation}
\sigma=\langle\tilde{\sigma}\rangle=\mathbb{C}<\tilde{\varepsilon}\rangle=\mathbb{C} \varepsilon
\end{equation}

In particular, by applying unit deformation conditions on “1,” 2 and 12 we obtain respectively the first, second and third column of C. So it is possible to do this in three separate load processes and consider the overall contribution as we are in linear elasticity under the hypothesis of small deformations for which the superimposition of the effects holds. The numerical simulation campaign is solved through the finite element method using triangular elements with three nodes. This is implemented through Mathematica software using the AceFEM finite element tool. The TEST [] procedure, written ad hoc, is used, which allows the three displacement conditions to be applied in sequence and thus obtain the stiffness matrix:

\begin{equation}
\mathbb{C}^{(o m)}=\left(\begin{array}{ccc}
756.78 & 80.90 & 0.00 \\
80.90 & 756.78 & 0.00 \\
0.00 & 0.00 & 40.87
\end{array}\right)
\end{equation}

From this matrix it is possible first of all to verify that the lower and upper limits (Reuss and Voigt) are satisfied. Furthermore, an isotropic behavior is observed and how the high pore fraction and therefore the large presence of low-stiffness material lowers the stiffness of the bone matrix. Further information on the automation and convergence of this result can be found in the convergence analysis section §8.4.

Although this first homogenization result may seem sufficient to consider bone matrix and pores as a single entity with such material properties, it is good to check what happens if the properties of the single component or its spatial arrangement vary. The following steps attempt to investigate whether this happens and how it affects overall behavior.

Microstructure variability

The load on the spine falls mostly on the vertebral body and the intervertebral discs. The vertebral body is formed of high porosity trabecular bone with a highly dense outer shell. However, the outer shell is very thin and is formed by a thickening of very compact trabecular bone, histologically different from the cortical bone. Several finite element analyzes have estimated that the outer shell contributes less than 15% to the overall load capacity of the vertebra [4, 11].

Furthermore, different trabecular bone samples showed different mechanical properties depending on the region where they are collected [12, 13]. This is a sign of the response of the bone which, being a living tissue and in continuous remodeling, is variable and adapts to the environment. The increase in compressive strength properties in the central part of the bone is nothing other than the response to the greater vertical load transmitted by the pulpy part of the vertebral disc compared to the adjacent peripheral region of fibrous rings.

Based on these considerations, further analyzes were conducted by varying the geometry of the RVE and observing how the resulting material response varies. The variations considered follow from the idealization of some physiological phenomena that lead to bone remodeling induced by the mechanical environment in which it is immersed.

Symmetry variations

The trabeculae therefore respond to the vertical load by positioning themselves mainly along the lines of force of the load and increase their density and thickness. An RVE with variable pores was considered to investigate the variation in material properties. A simple but accurate way to analyze this behavior is to consider elliptical pores. By introducing elliptical pores, the response of the bone to the vertical load can be interpreted as a reduction in the axial ratio and therefore a narrowing of the ellipse along the x direction.

The simulation analysis conducted show different changes in the stiffness matrix and therefore in the mechanical properties. Starting from a very low axial ratio, therefore a very flattened ellipse, there is a large increase in parameter C22 and a slight increase also in parameter C11 compared to the reference case (§4.3). The material symmetry and isotropy of the stiffness matrix is completely lost, a sign that the stiffness response along the direction of the larger half-diameter is greater than that in the direction of the smaller half-diameter. This reflects how much the bone produces with its remodeling response.

With an axial ratio of AxR = 0.1 the C22 coefficient reaches 5.376 GPa against a C11 = 1.462 GPa. This difference is reduced by decreasing the axial ratio. The two coefficients converge in the case in which the axial ratio tends to 1 or in the case in which C11 = C22 = C (om) 11. The Reuss and Voigt limits are verified in all cases. The data are present in fig. 3 and fig. 4. In decreasing the axial ratio, the reduction of the stiffness coefficients follows a linear trend. The C22 coefficient is significantly affected by this effect showing a much higher rate of decrease than the C11 coefficient which instead follows an almost constant trend.

It remains evident that as the axial ratio decreases, and therefore also the volume fraction, both the Reuss and Voigt limits and the numerical homogenization show the higher coefficient tending to the C (b) 22 coefficient of the bone matrix alone, while remaining lower at its 80%.

The same reasoning can be made in the opposite case in which diseases, such as osteoporosis, lead to a reduction in bone mass. This, revisited with a view to an increase in the axial ratio, shows how trabecular bone tends to reduce its structural properties even in areas where previously bone remodeling may have previously brought about a structural improvement (fewer pores and greater resistance). Morphological changes including the thinning of the trabeculae, an increase in the inter-trabecular space and the loss of connectivity can lead to a drastic decrease in resistance up to the risk of vertebral fracture [4].

Extension of the survey area

A further investigation was conducted by analyzing if and how the mechanical response varies by varying the window of interest. Having idealized a symmetrical RVE with symmetrical pores arranged on a regular grid it is clear that simply increasing the side would only include more pores. The volume fraction will remain the same as will the mechanical response. Even if the sum of the stresses increases, they are averaged over a larger area and their ratio, therefore the homogenization matrix, remains constant. The study then moves to a differently idealized RVE. The cases of interest are shown in fig. 5. The analysis is divided into two cases:

• Response analysis mechanics with a causal arrangement of pores according to a precise distribution identified by Doktor et al. [6].

• Mechanical response analysis with a hypothesized pore distribution starting from a real case, patient specific [14].

In analyzing the pores arranged according to a known distribution, different inclusions of fixed radius r = {0.1, 0.15, 0.2, 0.03, 0.35} are considered. These pores are randomly arranged by first considering 100, fig. 5 (a), and then 300, fig. 5 (b). For the computational implementation more information in §8.3.4.

We obtain coefficients of the stiffness matrix C11, C22 higher than the case §4.3. This is partly due to a lower pore fraction, i.e. a greater presence of bone matrix which contributes to increasing overall stiffness. It is interesting to note that the contribution of the volumetric fraction of pores is reflected in a greater stiffness along the main directions even if the structure does not seem to have a preferential direction for bone growth as in §4.3.

To reinforce these conclusions are also the results of the previous analysis (§5.1). Comparing them with these data we can see how the dominant factor remains the pore fraction but also as if the quantity of bone matrix disposed along the load direction increases, the stiffness in that direction increases.

In increasing the region of interest, 300 pores are considered and the side of the RVE is increased. This allows to include a greater number of pores and their distribution is such as to slightly reduce the coefficients of the stiffness matrix. The exact symmetry of the stiffness matrix is ​​lost but relatively stable coefficients C11 and C22 remain. The effect is less evident in the case of a greater number of inclusions. The loss of the isotropy properties of the matrix is ​​very evident in the coefficients that one would expect to be zero but even if very small, they are numbers other than zero. These latest results differ from homogenization in the ideal case and show overall greater stiffness attributable to the increase in bone mass.

These numerical results are influenced by the spatial distribution of the pores and different arrangements can lead to slightly different stiffness matrices while remaining below the limits of the isotropy loss observable in the case §5.1.

Specific case analysis

A further analysis of interest, similar to the previous case, considers a pore arrangement identified starting directly from the mentioned Micro CT image set. Taking a reference slice and selecting a specific RVE, it is possible to extract a specific bone structure fig. 5 (c) and (d). Manually designed circular pores were considered such as to be tangent to the bone matrix and to satisfy the pore fraction compatible with the data available in the literature (about 60%).

The first evident result of these tests is the increase in stiffness clearly due to the lower pore fraction considered, or rather to the increase in bone mass. The coefficients C11 and C22, although numerically different, remain very similar but larger than the previous cases but lower than the borderline cases found with the reduction of the axial ratio.

Material parameters

For the first homogenization study, average material parameters were considered among those identified in the literature. However, it is also known that the elastic modulus of bone structures can vary in relation to age or particular pathologies. Several results show that young people have higher stiffness modules than older people [2].

Both the contribution of the elastic modulus of the bone mass, Eb, and of the pores, Em. In particular, the increase in the elastic modulus of the bone mass gives more and more stiffness and a constant increase of the C11 parameters. As Eb varies from less than 2 GPa to slightly more than 6 GPa, the C11 coefficient goes from just under 300 to 1000. In the middle we find the case dealt with in the first study (4.3). A polynomial fitting satisfied by the relation:

\begin{equation}
\mathbb{C}_{11}=\mathbb{C}_{11}\left(E_{b}\right)=0.184+0.151 \cdot E_{b}
\end{equation}

Conversely, by varying the pore module, Em, no major changes are observed. Varying from 0.15 MPa up to 1000 MPa, a variation of less than 1% is observed. This reflects the fact that the greatest contribution of stiffness is given by the bone matrix. At the same time it tells us that, on a numerical level, it is also possible to use a module slightly higher than hypothesized in order to avoid too small numbers that can give computational problems going beyond the numerical precision of the calculator.

Conclusions

In light of the results obtained, it is clear that the predominant element in influencing the overall stiffness is the fraction of pores. In carrying out the homogenization, the fraction of pores can be taken into account and a reliable numerical result can be obtained.

However, it is good to consider that there is also a strong influence of the axial relationship in the pores. Since the bone is a living and continuously remodeling structure, it is also good to take into account the preferential directions of load along which a reduction in the axial ratio could be observed. Neglecting this effect and considering the trabecular bone as isotropic could lead to a homogenization matrix with significant numerical errors.

A strong weight is also given by the stiffness modules of the individual constituents but the greatest weight is given by the bone matrix which is responsible for giving stiffness to the overall structure. It is possible to avoid a specific case study and consider parameters and the error remains within acceptable limits. In some cases, however, this could underestimate the stiffness.

Possible future developments may be to analyze the influence of additional parameters such as trabecular thickness or connectivity. In addition, one could move on to a more advanced model of trabecular bone by considering the inside of the pores no longer as a highly elastic material but filled with fluid with the possibility of flowing and considering a poroelastic model.

Reference

1. Monesi V. Istologia. Piccin; 2012.
2. Cowin S. Bone Mechanics Handbook. CRC Press; 2001.
3. Nazarian vSDZDMRSBD A. Bone volume fraction explains the variation in strength and stiffness of cancellous bone affected bymetastatic cancer and osteoporosis. Calcified tissue international 007;83:368– 379.
4. Ferguson SJ, Steffen T. Biomechanics of the aging spine. European spine journal 2003;12:S97–S103.
5. Hamed E, Jasiuk I, Yoo A. Multi-scale modelling of elastic moduli of trabecular bone. Soc Interface
   2012;9.
6. Doktor T, Valach J, Kytyr D, Jiroušek O. Pore Size Distribution of Human Trabecular Bone – Comparison of IntrusionMeasurements with Image Analysis p. 115–118.
7. Aboudi J, Arnold SM, Bednarcyk BA. Chapter 3 Fundamentals of the Mechanics of Multiphase Materials. In: Aboudi J, Arnold SM, Bednarcyk BA, editors. Micromechanics of Composite Materials Oxford: Butterworth-Heinemann; 2013.p. 87–145. https://www.sciencedirect.com/science/ article/pii/B9780123970350000033.
8. DalstraM, Huiskes R, Odgaard A, Van Erning L. Mechanical and textural properties of pelvic trabecular
   bone. Journal of biomechanics 1993;26:523–35. https://doi.org/10.1016/0021-9290(93)90014-6.
9. Wirtz DC, Schiffers N, Pandorf T, Radermacher K, Weichert D, Forst R. Critical evaluation of known bone material properties to realize anisotropic FEsimulation of the proximal femur. Journal of
   biomechanics 2000;33:1325–30. https://doi.org/10.1016/s0021-9290(00)00069-5.
10. Jansen LE, Birch NP, Schiffman JD, Crosby AJ, Peyton SR. Mechanics of intact bone marrow. Journal
    of the mechanical behavior of biomedical materials 2015;50:299–307. https://doi.org/10.1016/j.
    jmbbm.2015.06.023.
11. Silva MJ, Keaveny TM, HayesWC. Load sharing between the shell and centrum in the lumbar vertebral body. Spine 1997;22,2:140–50.
12. Keller TS, Hansson TH, Abram AC, Spengler DM,PanjabiM. Regional variations in the compressive properties of lumbar vertebral trabeculae. Effects of disc degeneration. Spine 1989;14:1012–1019.
13. Keller TS, Ziv ME I, Spengler D. Interdependence of lumbar disc and subdiscal bone properties: a report of the normal and degenerated spine. Journal of spinal disorders 1993;6:106–113.
14. Beller G, Burkhart M, Felsenberg D, Gowin W, Hege HC, Koller B, et al. Vertebral Body Data Set ESA29-99-L3 2005;http://bone3d.zib.de/data/2005/ESA29-99-L3/.
15. Mastrofini A. Image processing volumetrico su osso trabecolare. website 2021 ; https://alessandromastrofini.it/image-volumetric-bone/.
16. Cowin S. Integrated Bone Tissue Physiology. Bone Mechanics Handbook 2001;1:4–5.
17. Mathworks. Find circles using circular Hough transform 2012;https://it.mathworks.com/help/images/
    ref/imfindcircles.html/.

Download the report

Omogeneizzazione numerica di materiali microstrutturati

Download Now!

GitHub repository