Structural analysis of bone by segmentation and finite element analysis in patients with osteoporosis

In bioengineering, mechanical properties related to the stress and stiffness of the materials can be determined by evaluating the bone density obtained from medical images. Diseases such as osteoporosis involve the decrease in bone mass and average bone density, inducing a high risk of fracture and the decrease in mechanical properties. This work proposes the use of numerical methods and medical imaging to evaluate biomarkers related to bone stress and strain affected by osteoporosis. The digital twin of bones is made from the segmentation of the bones associated with the femoral group diagnosed with osteoporosis. The material model for the bone tissue is defined as a functionally graded material with elastic properties obtained from computed tomography numbers. The results show that the proposed method allows monitoring of the disease and the evaluation of the risk of fracture, considering the functional mechanical load. In this sense, the use of radiomics takes advantage of the potential of personalized medicine to improve the lives of patients.


Introduction
Osteoporosis is defined as a systemic skeletal disease characterized by low bone density and deterioration of the microarchitecture of bone tissue with the consequent increase in bone fragility [1]. The measurement of bone mineral density (BMD) is commonly used as a clinical indicator of osteoporosis and fracture risk [2]. The decrease in BMD causes low skeletal strength, modifying the natural condition of bone properties in a healthy state, causing a high index of brittleness that in turn increases the susceptibility to fractures [3].
According to the World Health Organization (WHO), three million people suffer hip fractures each year, of which 200,000 die [3], and a third is registered in Latin America and Asia. It is believed that the structure or spatial arrangement of bone at the macroscopic and microscopic levels provides additional and independent information and may help to better predict the probability of fracture, in addition to evaluating the response to pharmacological intervention [4].
The main objective of image processing is the demarcation of objects that appear in digital images; this process is called segmentation, and a good approximation can often be achieved through thresholding. Generally speaking, this involves separating radiolucent areas from radiopaque ones and thereby identifying hypodense over hyperdense segments [5]; the attenuation coefficient of bone tissue is much higher than that of surrounding soft tissues, resulting in well-contrasted edges [6]. However, 2 apart from the surface description, more information can be obtained from the internal composition of the tissues using the Hounsfield (HU) scale [7,8].
On the other hand, numerical modeling allows approximating solutions to a large number of physical phenomena characterized by various types of complexities, such as geometric, material, or boundary conditions [9,10]. Finite element analysis (FEA) provides a better description of the stress and strain fields in a bone structure, this being a key factor to understand the functional behavior of the bone [11,12]. A digital twin of the bone tissue can give insight into the response of the tissue to different loading conditions over time, as the disease progress.
In this work, we investigate the feasibility of digital twins to assess the condition of bone tissues affected by osteoporosis. We perform bone segmentation and map the elastic properties obtained directly from the HU units. Then, the mechanical response is evaluated to obtain stresses and strains, and to estimate the risk of fracture.

Materials and methods
The three-dimensional (3D) model is developed from the computed axial tomography (CAT) of a patient between 40 and 50 years old, male, and of Colombian nationality. The images were provided by the High Technology Medical Center (CATME). The methodology is composed of four stages.

Model definition using segmentation
For the construction of the 3D computer-aided design (CAD) model, we used the 3D Slicer segmentation software, open-source software that allows the visualization and processing of medical images. This software has a variety of algorithms for the application of different segmentation techniques [13]. In addition, it offers tools for the refinement of the surfaces of the models.
Bone structures appear in two forms of tissues: cortical and trabecular [14]. To facilitate the visualization of the tissue and to identify the cortical and trabecular structure of the bone, we cut the CAT scan using the crop volume tools, Figure 1. We can separate unwanted tissue for processing, as seen in Figure 2. The manual segmentation process is performed using Grow from Seeds approach and Region Growing, which is a semi-automatic reconstruction method based on the evolution of seeds located by the user, to generate regions with properties similar to these seeds.
The most difficult part of the segmentation process occurs at the interface of the bone with other tissues, resulting in poor reconstruction of the surfaces. 3D Slicer offers options to smooth and create a uniform surface on the outside of the bone model. The three-dimensional model obtained in the 3D Slicer software is a digital representation of the geometry of the femur that does not have any associated information on the mechanical properties of the bone. The information is exported in stereolithography (STL) format for later meshing in Ansys.

Finite element model
We consider a linear elastic model with a standard Galerkin formulation. For the material model, we consider the bone as a functionally graded material (FGM), with elastic properties dependent on the position. This allows us to consider the anisotropic behavior of the bone tissue. The process for finite element analysis begins with the creation of a mesh from the CAD model in ANSYS software. Considering the complexity and precision that must be had for the calibration of the material properties, quadratic tetrahedral elements were selected for the mesh, which allow a good estimate of the apparent elastic modulus and provide better precision in the evaluation of the mechanical properties [15].

Calibration of apparent density and elastic modulus
To define the mechanical properties of the bone using information extracted from the CAT scan, the Bonemat v3.2 software was used, which allows assigning to each finite element (FE) the elastic properties of the bone from the information available in the diagnostic medical image. The apparent density ρ !"" is linearly related to the Hounsfield unit HU obtained from the CAT in Equation (1).
The HU unit refers to the quantitative scale used in computerized axial tomography to describe the different levels of radio-density in human tissues, a and b are the linear constants to calculate the apparent density of bone [16]. In this work, the values of the coefficients for calibration a = -0.021075 and b = 0.000786 were used, with the computed tomography (CT) number of air -1000 HU and water 0 HU. The theoretical and empirical relationships used for the definition of the material model considered a radiographic density of 1840 HU, for an elastic modulus value of 22000 MPa and apparent density of 1.73 g/cm 3 [17]. The elastic modulus E is obtained from the relation in Equation (2).
where a, b, c are calibration constants. For the process carried out in the Bonemat software, information from the DICOM file used in the segmentation is mapped to the elements in the mesh of the bone model.

Boundary conditions
We apply the displacement and load conditions necessary to perform the static analysis of the bone structure, where the external forces acting on the femur correspond to the state of equilibrium during standing in unipolar support [18], as shown in Table 1 for the force applied on the femoral head, and Table 2 for the force applied on the greater trochanter.  This study was carried out for a person with a mass of 100 kg, from this mass the forces acting on the femur were calculated. Moreover, the study was carried out for the bone assuming conditions of a homogeneous material in the sections of trabecular bone and cortical bone, assigning constant average values of elastic modulus and density shown in Table.   Table 3. Material properties [19].

Results
The finite element model for the femur was created with 61201 solid elements of the quadratic tetrahedron type. The assignment process of material properties for the bone with osteoporosis produced a total of 573 materials, with graduated properties of elastic modulus and densities, distributed in the elements of the model. Equivalent von Mises stress and principal stress distributions were obtained, which can be compared with the elastic limit of the material. Considering the behavior of the bone as that of a brittle material [20], the results of the principal stresses are taken into account. Table 4 shows the maximum von Mises, principal stress, equivalent strain, and displacement for a healthy bone and with osteoporosis. Notice that the maximum stress for bone with osteoporosis is greater than for bone with a healthy condition. This behavior occurs because in bone with osteoporosis the stress increases with the decrease in density and bone loss. In this way, the affected bone tends to reach fracture more easily in the area of the femoral shaft, which is where the maximum point of stress is presented.
The maximum strains are located near the femoral head for both, the healthy bone, Figure 3(a), and with osteoporosis, Figure 3(b). The maximum value of stress is presented in the diaphysis, and the stress is distributed throughout the cortical tissue until reaching the lower part of the proximal epiphysis (neck of the femur), Figure 3. The maximum stress of bone with osteoporosis, Figure 3(c), shows an increase of 24.71% with respect to healthy bone, Figure 3(d), and both are located in the same region.

Conclusions
The segmentation process using 3D Slicer provides computer-aided design models of bone structures, separated from other tissues. The gradual assignment of elastic properties from the computed tomography number allows evaluating the mechanical behavior of bone structures affected by osteoporosis, unlike homogenized models. Osteoporotic bone exhibits less stiffness and greater stresses compared to unaffected bone, due to the loss of both cortical and trabecular bone, evidenced by the loss of density. The proposed methodology using a digital twin of the bone tissue allows monitoring the disease and assessing the risk of fracture with greater precision considering functional mechanical loads. The use of radiomics leverage personalized medicine, focusing on the profile of the patient instead of the disease.