Simulating Bone Atrophy and Its Effects on the Structure and Stability of the Trabecular Bone

According to the world health organization (WHO) osteoporosis is considered to belong to the ten most important diseases worldwide. It is defined as a skeletal disorder characterized by compromised bone strength predisposing to an increased risk of fracture (NIH, 2000). Osteoporosis is a metabolic bone disorder in which bones become brittle and prone to fracture. The underlying reason for the occurrence of osteoporosis is that the two processes being responsible for bone remodelling turn out of balance. Bone tissue is continuously resorbed by special bone cells, the so-called osteoclasts. On the other hand, the formation of bone also takes place by the action of other bone cells, namely the osteoblasts. These cells form new bone. The bone formation is, however, not uniform. It is rather controlled by an external mechanical stimulus such that more bone material is produced at those sites where the local stress is larger. This leads to an adaptation of the inner bone structure (trabecular bone) to externally acting forces on the bone (Mullender & Huiskes, 1995). Thereby, a minimal-weight structure, that is adapted to its applied stresses, is formed as it was already correctly conjectured by Julius Wolff as early as in 1892 (Wolff, 1892). In a healthy bone there exists an equilibrium between bone formation and bone resorption, whereas in an osteoporotic bone more bone resorption than formation takes place, which leads to a rarefied network of the trabecular bone. Besides this disease induced effect the rarefication of the inner bone structure can also have other causes like living under zero gravity conditions, immobility or age-related atrophy as the most common example. Advances in modern imaging modalities like high resolution magnetic resonance (HRMR) or micro computed tomography (μCT) imaging have led to great improvements of the image quality especially in terms of spatial resolution which now allows for a proper threedimensional visualisation of the trabecular network. Having these highly resolved data at


Introduction
According to the world health organization (WHO) osteoporosis is considered to belong to the ten most important diseases worldwide. It is defined as a skeletal disorder characterized by compromised bone strength predisposing to an increased risk of fracture (NIH, 2000). Osteoporosis is a metabolic bone disorder in which bones become brittle and prone to fracture. The underlying reason for the occurrence of osteoporosis is that the two processes being responsible for bone remodelling turn out of balance. Bone tissue is continuously resorbed by special bone cells, the so-called osteoclasts. On the other hand, the formation of bone also takes place by the action of other bone cells, namely the osteoblasts. These cells form new bone. The bone formation is, however, not uniform. It is rather controlled by an external mechanical stimulus such that more bone material is produced at those sites where the local stress is larger. This leads to an adaptation of the inner bone structure (trabecular bone) to externally acting forces on the bone (Mullender & Huiskes, 1995). Thereby, a minimal-weight structure, that is adapted to its applied stresses, is formed as it was already correctly conjectured by Julius Wolff as early as in 1892 (Wolff, 1892). In a healthy bone there exists an equilibrium between bone formation and bone resorption, whereas in an osteoporotic bone more bone resorption than formation takes place, which leads to a rarefied network of the trabecular bone. Besides this disease induced effect the rarefication of the inner bone structure can also have other causes like living under zero gravity conditions, immobility or age-related atrophy as the most common example. Advances in modern imaging modalities like high resolution magnetic resonance (HRMR) or micro computed tomography (µCT) imaging have led to great improvements of the image quality especially in terms of spatial resolution which now allows for a proper threedimensional visualisation of the trabecular network. Having these highly resolved data at 696 hand it now becomes possible to perform a differential analysis of both the structural and the mechanical properties of the trabecular bone. We determine the local structure of the trabecular network by calculating isotropic (α) and anisotropic scaling indices (α x ,α y ,α z ) (Räth et al., 2008). These measures have been proven to be able to discriminate rod-from sheet-like structures and to quantify the alignment of structures with respect to a preferential direction as e.g. given by the direction of the external force. Another class of texture measures is given by Minkowski Functionals (Michielsen & De Raedt 2001), which -as the scaling indices -also incorporate correlation functions of higher orders and supply global morphological information about structures under study. The calculation of the local mechanical properties, i.e. the load distribution inside bodies with high porosity and complex architecture is enabled by the use of the finite element method (FEM), which has become a standard tool in bone research. Bone modelling and remodelling processes can be simulated by describing the action of the osteoblasts and osteoclasts by means of rate equations. The equations are nonlinear partial differential equations, which can thus only be solved numerically (Huiskes et al, 2000). As it is mostly the case in nonlinear systems the solutions depend very sensitively on the choice of parameters of the model, with which e.g. the onset of bone formation is controlled. In this chapter we propose methods to simulate the effects of bone modelling on the structure and stability of the trabecular bone. We gradually deteriorate the trabecular structure by using some concept of cellular automata (Wolfram 1983, Wolfram 1984 to solve the rate equations numerically. Specifically, both the three space coordinates as well as the time is discretised. The bone resorption, which is described as a continuous process by the partial differential equations, is thus transformed to a consecutive removal of bone surface voxel, which is discrete in space and time. Having simulated bone modelling we then assess the effects of changes in the bone structure on both the structural and mechanical properties of the specimen. We do not consider the details of the physiological mechanisms of the bone adaptation to mechanical loading and do not describe metabolic activity of the biological cells. Rather, we assume independent action of osteoclasts (bone resorbing cells) and osteoblats (bone forming cells), which is typical for bone modelling process leading to global morphological and topological changes. Among the large variety of bone modelling scenarios we concentrate on bone atrophy due to erosion of the bone surface by osteoblasts resorption activity and decrease of osteoblast formation activity, which is typical for the process of normal aging (often called primary osteoporosis). This type of bone modelling can be simulated by random resorption of bone voxels at the interface between trabecular bone and bone marrow. We develop three numerical models for simulation of bone loss, which correspond to the different assumptions about age-related atrophy in male and female bones: thinning of the trabecular bone structure with and without preserving of topological connectivity and with preferential loss of aligned rod-like trabecular elements, which were previously identified by a scaling index analysis.

The data set
For our study we chose 17 specimens of young (50 < age < 70) patients with high maximum compressive strength (MCS > 70 Newton (N)) out of a data set (Räth et al., 2008) of 151 cylindrical specimens with a diameter of 8 mm and a length of 14 mm, which were harvested from 73 thoracic and 78 lumbar human vertebrae. The resulting µCT grey-value images with isotropic spatial resolution of 26 µm were segmented using a low-pass filter by convolving the image with a Gaussian kernel with standard deviation 0.8 and support of 1 to remove noise and a fixed global threshold equal to 22% of the maximal grey value to extract the mineralised bone phase (Hildebrand et al 1999). After µCT scanning, the bone samples were cut to the length of 12 mm and tested by applying uniaxial mechanical compressive load using a servo-hydraulic machine (MTS 858 mini Bionix II, MTS Eden Prairie, USA) with a load cell of 1.5 kN. Maximum compressive strength (MCS) was determined in biomechanical experiment as the first local maximum of the forcedisplacement curve and used in correlation analysis as a golden standard for characterisation of bone strength (Eckstein et al., 2007). Main structural characteristics of the bone specimens are given in Table 1. The binarized 3D µCT images were used as a starting structure for the numerical simulations of the bone resorption process. According to integral geometry an n-dimensional body can be completely characterized by n+1 functionals, which evaluate both size and shape of the object. In a three-dimensional space the four functionals are represented by the volume (MF 1 ), surface area (MF 2 ), integral mean curvature (MF 3 ) and integral Gaussian curvature (MF 4 ). Minkowski Functionals are derived from the theory of convex sets and expressed as volume integral for MF 1 and surface integrals over boundary S of the excursion set I with the principal radii of curvature R 1 and R 2 for other functionals. The first two functionals MF 1 and MF 2 describe morphology of the structure and coincide with morphometrical parameters bone volume and surface fractions (MF 1 = BV/TV and MF 2 = BS/TV). The fourth integral, also known as Euler characteristic χ, characterises the topological connectivity of structures and can be expressed in terms of Betti numbers β 0 (number of connected components), β 1 (number of tunnels), β 2 (number of cavities): In the case of binary images we have exactly four global characteristics, which can be used as texture measures for assessment of bone strength (Monetti et al., 2011).

Isotropic and anisotropic scaling indices
We calculate isotropic and anisotropic scaling indices (SIM) (Räth & Morfill 1997;Monetti et al. 2004;Müller et al. 2006;Räth et al. 2008) as measures to characterize the complex trabecular network and its degree of alignment relative to a preferential direction. Generally, scaling indices represent one way to estimate the local scaling properties of a ndimensional point set. Considering binarised three-dimensional µCT-images, a suitable representation of the image information as a set of three-dimensional points is given by for the isotropic case. The calculation of the scaling indices depends on the parameters n and r. The exponent n controls the shape of the weighting function. With increasing n the weighting function becomes more and more step-like. In all our studies we found that n has only little influence on the results. For the following calculations we fixed n to n=2, which emphasises the connection to Gaussian smoothing or kernel functions. The scale parameter r controls to which scales of the image structure the scaling indices are sensitive. From previous studies we found that r = 8 (in voxel units) is a well-suited choice. By means of the scaling indices α one can then discriminate between voxels belonging to rod-like (α ≅ 1) and sheet-like (α ≅ 2) structural elements of the trabecular network, which is one of most important discriminating feature to discern healthy from osteoporotic bone structure. It is straightforward to implement anisotropies in the calculation of the scaling indices by introducing an ellipsoidal distance measure with the eigenvalues λ x , λ y and λ z : In this study we calculate anisotropic scaling indices α z with respect to the direction of the external force acting in the bone, for which we chose the ratio of 5:1 for the eigenvalues, i.e. λ z = 5λ x = 5λ y , for setting the degree of anisotropy along the z-direction. This setting was found to be well suited especially for the detection of aligned cylindrical structures (i.e. the trabeculae) with typical mean breadth and length.

Finite element models
For assessing biomechanical strength of trabecular network we apply Finite Element Method (FEM) with linear elastic assumptions (Rietbergen et al., 1995;Rietbergen et al., 1999;Sidorenko et al. 2009). For the bone mineral material we assume isotropic properties with Young's modulus Y = 10 GPa and Poisson's ratio ν = 0.3 and describe relationship between stress σ and strain ε components by generalized Hooke's law (constitutive equations), which states that stress is proportional to the strain up to the elastic limit We use Dirichlet boundary conditions with prescribed on the top surface constant strain ε 0 = 1% to simulate uniaxial loading in the natural direction (in this work denoted as z axis) applied in biomechanical experiment. We generate a finite element mesh by direct converting image voxels that belong to the bone tissue into equally sized and oriented eight node brick elements. An exact number of nodes depend on the structure characteristics and varies for our data set from 0.7⋅10 6 for weak bones, dominated by rod-like trabecular elements, up to 2⋅10 6 for strong bones with a lot of plate-like formations. From discrete nodal displacements obtained from FEM strain and stress components can be recovered at any point of the structure. For correlation analysis with respect to the MCS we use apparent total reaction force F r at the top face A t which is recalculated from the normal stress component in direction of applied mechanical load σ t zz :

Bone resorption models
We develop three numerical models for bone resorption, which correspond to the different assumptions about bone atrophy in male and female bones. Several investigations (Aaron et al., 1987;Khosla et al., 2006) demonstrated sex difference in trabecular bone aging. Although the decrease with age in trabecular bone volume is common for both sexes, male cancellous bone exposes uniform thinning, whereas female trabecular network suffers from loss of connectivity and entire rod-like structural elements (Barger-Lux et al., 2002). In all three models we simulate the random resorption of bone material by removal surface voxels according to the given value of relative bone loss volume ∆BV/TV, but with different topological features. In order to study the role of topological connectivity for bone strength in the first resorption model (Model I), a bone voxel is removed only if does not change topology of the system, i.e. no new cavities or tunnels are created. In terms of global topological characteristics it means that the fourth Minkowski Functional, which coincides with Euler characteristic χ is conserved, We compare this model with two other models, which both do not preserve the connectivity (χ ≠ const). In Model II there isn't any limitations or conditions on the removal of bone surface voxels. In Model III, however, we simulate the preferential destruction of rod-like trabecular elements by only removing surface voxels with local topological anisotropic scaling index α z < 2. Typical changes in trabecular bone structure are demonstrated in Fig. 1. Three regions in circles show the differences in structure due to different resorption models: preserving of connectivity and destroying of rod-like elements at the same place in the trabecular network.

Results
As a main characteristic of the statistical analysis we use Pearson's correlation coefficient r MCS (Tables 2 -5)  Plots for MF 4 = χ (Fig. 2, last row) proof that in first resorption model (left column) the connectivity is preserved (χ = const) and in the two other models the porosity of the structure increases (∆χ < 0) during the resorption process. In both models the increase of . In the third model with preferential resorption of rod-like trabecular elements (with χ ≠ const and α z < 2) the growth of number of tunnels (β 1 ) is compensated by an increase of number of separate parts (β 0 ) and slows down negative increase of χ. For both resorption models without conservation of connectivity for structure with large relative bone loss (∆BV/TV > 30%) the correlation coefficient of MF 4 = χ with MCS becomes considerably higher than that for original structure (Tables 2). www.intechopen.com  Table 2. Pearson's correlation coefficient of the four MF with respect to MCS for the three models of bone loss and eight rareficiation steps ranging from 0% to 50% removal of the initial bone volume.
For isotropic SIM (Fig. 3a) there is almost no difference in the P(α) spectrum for the different resorption models observed, which leads to only a small decrease in the correlation coefficient (last row in Fig. 3a and Table 3a). For anisotropic SIM (Fig. 3b) there is an obvious difference in P(α z ) already at ∆BV/TV = 10% and a more considerable decrease of r MCS (last row on Fig.3b and Table 3b) in the case of preferential resorption of rod-like trabecular elements (third model with χ ≠ const and α z < 2).
SIM α χ=const χ≠const χ≠const, α z <2 Fig. 3a. Probability distribution function of isotropic scaling index P(α) for 17 bone specimens with original structure (first row) and with different value of bone resorption (30%: second row, 50% third row). Last row: mean value of α spectrum (black curves and left axis) and correlation coefficient of mean value with MCS (red curve and right axis) for three models of bone loss (from left to right).

www.intechopen.com
Simulating Bone Atrophy and Its Effects on the Structure and Stability of the Trabecular Bone 705 SIM α z χ=const χ≠const χ≠const, α z <2 Fig. 3b. Probability distribution function of anisotropic scaling index P(α z ) for 17 bone specimens with original structure (first row) and with different value of bone resorption (30%: second row, 50% third row). Last row: mean value of α z spectrum (black curves and left axis) and correlation coefficient of mean value with MCS (red curve and right axis) for three models of bone loss (from left to right).  Table 3. Correlation coefficient of mean value of α spectrum (a) and α z spectrum (b) with MCS for three models of bone loss.
The mechanical strength of the resorbed trabecular structure as determined with FEM ( Fig.   4) depends almost linearly on the relative bone loss ∆BV/TV and shows small decrease in the correlation with MCS (Table 4). At high bone loss (∆BV/TV = 50%) the mechanical strength of all bone specimens was found to be larger in the case with conservation of connectivity (solid line in Fig. 5). In Table 5 we summarise effect of large bone loss (∆BV/TV = 40%) on different numerical texture measures. FEM and SIM demonstrate small drop in correlation with MCS of initial structure. These methods can be proposed for prediction of osteoporosis: relative strength and local topology do not change considerably under the process of random surface resorpton. In a contrast, after large bone resorption global MF 2 and especially MF 4 improve their correlation with MCS of the original structure. This effect can be used for the diagnostic of the current state of the bone structure.

Summary and conclusions
We proposed a method based on the ideas of cellular automata to simulate bone atrophy and applied it to a sample of 17 bone probes visualised with high resolution µCT imaging. Although our study is so far restricted to the simulation of bone loss and did not include any processes of bone formation, we could already gain some new and very interesting insights about the important factors determining the strength of bones. As expected we found that the initial structure determines the relative strength of the bone under random surface bone losses. Patients with stronger bones in young age have better prognosis for age-related bone atrophy. We found that the connectivity plays the most important role in determining the strength of the bone structure: among three resorption models the highest apparent reaction force was calculated for the resorption model which preserved the connectivity (χ = const). FEM, isotropic SIM, the first and second MF yielded stable values of the correlation coefficient r MCS under the random bone loss process for all numerical resorption models and can be recommended for prediction of bone strength in bone atrophy process.
The mean value of the anisotropic scaling indices α z demonstrated sensibility for preferential rod-like trabecular loss as simulated by third numerical resorption model (with χ ≠ const and α z < 2). For this model the scaling index approach shows a decrease of correlation coefficient r MCS already at 10% loss of mineral bone fracture.
For the two resorption models without conservation of the connectivity (χ ≠ const) the bone surface resorption significantly improves the correlation of the fourth MF with MCS (from r MCS = -0.12 for original structure up to the r MCS = -0.92 for bone loss ratio ∆BV/TV = 35%).
Such an effect suggests that the random surface resorption destroys thin and unimportant connections of the trabeculae and only the strong and thick trabecular elements are taken into account for correlation with MCS. The removal of bone voxels can thus be interpreted as a distillation of the essential skeleton of the trabecular structure, which is a much more sensitive tracer of the mechanical stability of the bone. In fact we found that fourth Minkowski Functional calculated for structure prepared by random surface resorption yields higher correlations with MCS than FEM-based measures, which are so far considered to yield the highest correlations with the mechanical properties of bone probes. Therefore the rarefication procedure as outlined in this study in combination with Minkowski Functionals may lead to a novel technique for the diagnosis of the trabecular bone quality and strength in the prediction of osteoporosis.

Acknowledgments
This study was in part supported by the Deutsche Forschungsgemeinschaft (DFG) under the grants MU 2288/2-2 and BA 4085/1-2. Osteoporosis is a public health issue worldwide. During the last few years, progress has been made concerning the knowledge of the pathophysiological mechanism of the disease. Sophisticated technologies have added important information in bone mineral density measurements and, additionally, geometrical and mechanical properties of bone. New bone indices have been developed from biochemical and hormonal measurements in order to investigate bone metabolism. Although it is clear that drugs are an essential element of the therapy, beyond medication there are other interventions in the management of the disease. Prevention of osteoporosis starts in young ages and continues during aging in order to prevent fractures associated with impaired quality of life, physical decline, mortality, and high cost for the health system. A number of different specialties are holding the scientific knowledge in osteoporosis. For this reason, we have collected papers from scientific departments all over the world for this book. The book includes up-to-date information about basics of bones, epidemiological data, diagnosis and assessment of osteoporosis, secondary osteoporosis, pediatric issues, prevention and treatment strategies, and research papers from osteoporotic fields.