Bone Stress-Strain State Evaluation Using CT Based FEM

Nowadays, the use of a digital prototype in numerical modeling is one of the main approaches to calculating the elements of an inhomogeneous structure under the influence of external forces. The article considers a finite element analysis method based on computed tomography data. The calculations used a three-dimensional isoparametric finite element of a continuous medium developed by the authors with a linear approximation, based on weighted integration of the local stiffness matrix. The purpose of this study is to describe a general algorithm for constructing a numerical model that allows static calculation of objects with a porous structure according to its computed tomography data. Numerical modeling was carried out using kinematic boundary conditions. To evaluate the results obtained, computational and postprocessor grids were introduced. The qualitative assessment of the modeling data was based on the normalized error. Three-point bending of bone specimens of the pig forelimbs was considered as a model problem. The numerical simulation results were compared with the data obtained from a physical experiment. The relative error ranged from 3 to 15%, and the crack location, determined by the physical experiment, corresponded to the area where the ultimate strength values were exceeded, determined by numerical modeling. The results obtained reflect not only the effectiveness of the proposed approach, but also the agreement with experimental data. This method turned out to be relatively non-resource-intensive and time-efficient.


INTRODUCTION
At the moment, numerical modeling is the main calculation method in various fields of scientific research. And the use of data from experimental samples' images is the most promising approach to describing the behavior of elements of inhomogeneous structures under external the influence (Mithun and Mohammad, 2014;Marwa et al., 2018). This approach allows us to consider problems with multi-connected porous or inhomogeneous materials (Kayumov, 1999;Alberich-Bayarri et al., 2008;Kharin et al., 2019). This is especially true in clinical practice since it can significantly affect the quality of treatment (Kichenko et al., 2008;Kichenko et al., 2011;Maquer et al., 2015). Stress-strain state is important in tasks to estimate bone remodeling, especially in the case of patient-based implant design (Maslov, 2017;Borovkov et al., 2018;Maslov, 2020). There are several image-based approaches to determining the mechanical parameters of the model. One of them is based on approximating the inhomogeneity by constructing a distribution using the mean intercept length method. With this approach, the physical relations are expressed in terms of the tensor of elastic constants and the fabric tensor (Sachenkov et al., 2016;Kayumov et al., 2018;Yaikova et al., 2019). The next approach is based on reducing the material anisotropy to orthotropy by calculating constants from numerical experiments (Ruess et al., 2012;Ridwan-Pramana et al., 2017;Chikova et al., 2018). Finally, the work presents a third approach based on taking into account the material properties by weighted integration of the global stiffness matrix of the finite element mesh.
Computed tomography is the main approach to reconstructing images of the investigated domain. Such a process implies the creation of its digital prototype, which is a three-dimensional array of numbers corresponding to the material permeability in each microelement (voxel) of the volume. These values can be interpreted according to Hounsfield's quantitative scale of X-ray density. Modern computed tomography (CT) allows getting detailed information not only about bone tissue state but also its interaction with implants (Stock et al., 2002;Jones et al., 2007;Vanlenthe et al., 2007). Using CT data the segmentation can be performed. It's a complicated process because it's important to restore not only geometry but also the inner structure of the organ (Supplementary Material). Wherein restored geometry should be meshed for simulation. The quality of mesh influences a lot on simulation quality. A common approach is to use stratified segmentation. E.g., separate geometry for cortical and cancellous bone tissue in whole bone (Marcián et al., 2018;Carniel et al., 2019). Then, in simulation, each geometry can be taken into account with different mechanical properties (Marcian et al., 2017;Hettich et al., 2018;Kuchumov and Selyaninov, 2020). But these actions are laborious and modern research focused on ways to facilitate the segmentation (Moriya et al., 2018). On other hand, the highest calculation accuracy is achieved by modeling each microvolume of a continuous medium with one finite element (Prez et al., 2014;Tveito et al., 2017). However, in this case, the computational intensity increases excessively. In this regard, it makes sense to increase the finite element size so that it contains a certain amount of microvolumes (Nadal et al., 2013;Marco et al., 2015;Giovannelli et al., 2017). In this case, each voxel can be considered as a neighborhood of the integration point of the local stiffness matrix. However, the question of the number of voxels in the element volume is still open. The same question is open for applicability of finite element method (FEM) in casual clinic practice (Viceconti et al., 2018).
The purpose of the study is to implement the method of static calculation of problems based on samples of a porous structure. A three-dimensional isoparametric linear element of a continuous medium, built on the basis of computed tomography data, was used as a finite element of the grid.

Computed Tomography Based Finite Element
Let's consider the region describing the volume of an object. Additionally, there is digital data of the object. In the research computed tomography data is understood as a digital prototype. The data consists of a number of microvolumes with an average value of the Hounsfield unit. Usually, such data illustrated by a greyscale three-dimensional image. Each microvolumes (voxel) has dimensions Δx, Δy, Δz. Usually, a voxel is equal to a cube and the size depends on computed tomography resolution. The mechanical behavior of the object will be investigated by the finite element methods. In this case, the volume is consists of a number of finite element. To simulate the inhomogeneity continuum three-dimensional eight-node finite element with bilinear interpolation was built. Using shape functions N n in the local coordinate system of the element (ξ n , η n , ζ n ) the displacement can be presented as follows: where u, v, w are components of a displacement vector in global coordinate system x, y, z, and shape function looks like N n (ξ, η, ζ) (1 + ξ · ξ n ) (1 + η·η n ) (1 + ζ·ζ n )/8. The equations can be presented in matrix from: where {ϑ E }-vector of the node displacements, [N]-matrix of the shape functions. Let's introduce strain and stress vectors: Using strain-displacement equations strain vector can be determined (Zienkiewicz and Zhu, 1987): Using Hooke's law stress vector can be found: where [D]-elasticity matrix, {r}-vector in a point where elasticity matrix set. In case of inhomogeneous media elasticity matrix can be rewritten in form: where [D B ]-elasticity matrix of pure material, ω{r}-scalar function which determines inhomogeneity by digital prototype data.
The ω{r} function can be calculated using Hounsfield unit (HU) from CT data. And then strain vector can be rewritten as follow: So stiffness matrix for an element can be denoted (Zienkiewicz and Zhu, 1987): where |J(ξ, η, ζ)| is the Jacobian determinant of coordinate transformation. When node displacements are found using differentiation matrix and elasticity matrix strain and stress vectors can be found. Using shape functions for every element stress function can be approximated: So average stress in each element can be found: To estimate computational error the energy norm should be calculated. For this purpose at every element stress error should be calculated: where {σ a n }-average stress vector at node n, {σ i n }-stress vector of node n of finite element i.
Then, for each element, the energy error is calculated as follows (Grassi et al., 2012;Giovannelli et al., 2017): Strain energy of the element can be calculated by equation: Finally, the energy error can be normalized according to the strain energy:Ẽ So, the described approach allows taking into account the inhomogeneity of media in terms of the stiffness matrix. Analyses of stress-strain state should be held in terms of average stress and strain. The quality of results can be judged by normalized energy error value. A more detailed description of the method and evaluation of element convergence depending on voxel number is given in the previous publications Gerasimov et al., 2021).

Meshing and Segmentation
Consider the global algorithm of the analysis of the bone. Input data is presented as CT, where each value of HU can be interpreted as optical density and even can be recalculated to elastic constants, and ultimate stress using equations (Rho et al., 1995): where coefficients a and b with the corresponding indices are determined from experiments or literature and can vary depending on the type of bone tissue or organ. Bone is a heterogeneous porous material with an irregular structure. The process of constructing a finite element mesh over the selected area implies segmentation of the initial CT data. In this case, each finite element is associated with a certain area of the bone. Such subregions are characterized by a different distribution of bone substance, which can include both the cortical plate and the spongy substance. Thus, the proposed method for constructing the finite element mesh makes it possible to take into account the material anisotropy at the level of modeling an individual finite element based on its CT data. Figure 1 sequentially shows subvolumes with different levels of porosity. On the left side is the area with a highest porosity (cancellous zone), on the right one is the part corresponding to a cortical plate. A popular approach is to use segmentation to restore the geometry of an organ. But researchers face difficulty when investigating organs with irregular geometry and complicated inner structure, e.g., femur, where there are cortical and cancellous regions. The proposed method allows taking into account such anisotropy properties indirectly. Based on the developed FE it is offered to use regular mesh for the whole CT space. Using a threshold filter for porosity the mesh can be clarified. Let's call such a mesh computational grid. Then kinematic or/and static boundary conditions can be applied. Then all local stiffness matrixes can be calculated using CT data. Assembling the global stiffness matrix is a casual task. After nodal displacements are found average stress and strain vector for each element can be found. Principal values, directions, and von Mises yield was calculated for average stress and strain in each element. To juxtapose stress-strain state and CT data additional mesh was build. Let's call such a mesh postprocessing grid. Such a grid can be built as a contour of binarized CT data. So, results from the computational grid can be interpolated to the post-processing grid.
General algorithm Input: CT data, boundary conditions Output: Post processing grid, nodal and element solutions 1: generate regular mesh for hole CT space with given element size 2: for each element in regular mesh 3: if bone tissue relative content < bone threshold 4: then delete element 4: else compute local stiffness matrix according to CT data Frontiers in Mechanical Engineering | www.frontiersin.org June 2021 | Volume 7 | Article 688474 4: End for 5: Save mesh as a computational grid 6: Assembling the global stiffness matrix, apply boundary conditions and solve 7: Calculate stress and strain in elements 8: Generate post-processing grid using given HU threshold with given tolerance 9: Interpolate solution to post-processing grid FEM software was implemented using C++ (g++ compiler). An Eight-node finite element with three degrees of freedom in each node was realized. The algorithm was parallelized using OpenMP technology because of the high labor costs of integrating the local stiffness matrix. The solution of problems was carried out on a computer with the following configuration: CPU-AMD Ryzen 7 1700 with 8 physical and 16 logical cores with a frequency of 3.7 GHz, RAM-G. Skill Aegis 16 GB DDR4 3000 16GISB K2 C16 with a frequency of 2933 MHz, motherboard-MSI B350M MORTAR. Post-processing was carried out using the ParaView software (Ahrens et al., 2005;Ayachit, 2015).

Bending Test
Experiments were performed on six 15-20 kg male Vietnamese swine. Animals were placed in individual cages under standard laboratory conditions with unlimited access to food and water and a 12 h day/night cycle. The protocol of the experiment, including anesthesia, surgery, postoperative care, testing, and euthanasia, was approved by the Animal Care Committee of Kazan State Medical University (protocol #5 of 20 May 2020). All experimental procedures were performed in accordance with the standards to minimize animal suffering and the size of experimental groups. Animals were included in experiments after a period of adaptation of at least 7 days. In experiments bones of animals after contusion spinal injury were used. The animals were sacrificed on the 8th week after a spinal cord injury (SCI). The euthanasia was performed under deep anesthesia by overdosing on an inhalational narcosis agent (Isoflurane) (Fadeev et al., 2020). The bones of the fore and hind limbs were extracted.
Bones were scanned by CT. Micro/nanofocal X-ray inspection system for CT and 2D inspections of Phoenix V | tome | X S240 was used for scanning. The system is equipped with two X-ray tubes: microfocus with a maximum accelerating voltage of 240 kV power of 320 W and nanofocus with a maximum accelerating voltage of 180 kV power of 15 W. For primary data processing and creating a volume (voxel) model of the sample based on x-ray images, the datos|x reconstruction software was used. The sample fixed in the holder was placed on the rotating table of the X-ray computed tomography chamber at the optimum distance from the X-ray source. The survey was carried out at an accelerating voltage of 90-100 kV and a current 140-150 mA. After that three-point bending tests were carried out.
CT data was about 800-900 Mb, voxel size in the direction of three Cartesian coordinate axes was equal to 0.2 mm. The data dimension was equal to 7523. Based on CT data, an initial regular finite element mesh was constructed. The computational grid generation was based on the removal of elements with a bone substance content of less than 5%. In this case, the initial CT data was pre-binarized according to the specified threshold value: values less than the threshold value were equated to zero, larger or equal-to one. Determination of the mechanical properties of bone tissue was based on the use of a similar binarization function. Young's modulus for pure bone tissue can be calculated by Eq. 2.2.2 and was equal to 30 GPa and Poisson's ratio equal to 0.3. In calculations three-point bending simulation. The top and the bottom surfaces of the distal and the proximal regions were fixed (marked by black lines in Figures  2A,C). Displacements in y negative direction were applied on the top surface of the medial region (marked by red arrows in Figures  2A,C) and were equal to 2, 3, and 4 mm, respectively. To verify the simulation results equivalent force was calculated. All steps of calculations from reading CT data to interpolate nodal and element solution to the post-processing grid took about 15 min.

RESULTS AND DISCUSSION
Two initial meshes with 10 and 5 mm element sizes were used. Normalized errors were analyzed for each result. Special attention was paid to the value of normalized errors in elements where maximum Von Mises stress appears. In Figure 2 distribution of normalized error (on computational grid) and Von Mises stress (on post processing grid) in forelimb is shown. It's important to understand correctly the meaning of the Von Mises stress in results. Averaging in element according to corresponding CT data change the understanding of the received value. In this case, we consider the element as a subvolume of the digital prototype with unknown stress spatial distribution but a known mean value. It means that falsely high and falsely low values can appear, that's why normalized errors should be checked in the element. For both meshes maximum of normalized error appears in border elements where relative bone tissue content is low. It can be Let's focus on the task with 4 mm displacement applied and consider the mesh with a 5 mm element size. The minimum of normalized error (20%) and maximum of Von Mises (500 MPa) stress was in the region near applied loading (see Figures 2C,D). For more detailed analyses principal stresses and directions were found, which are shown in Figure 3. To estimate equivalent force in displacements' applied region 2nd principal stresses were used (see Figures 3C,D). The integral of the stress vector in the region multiplied by the area will consider as a reaction to equivalent force. 1st and 3 days principal stresses describe the stress state of the bone (see Figures 3A,B,E,F). Corresponding streamlines form an orthogonal system. It should be noticed that in the region of Von Mises maximum first principal stresses reach a maximum value 400 MPa and 3 days principal-minimum −400 MPa. It explains the crack that appears in an experiment, cause the values exceed yield strength (Crenshaw et al., 1981;Imai, 2015).
To verify simulation results equivalent force was compared with force from three-point bending experiments. In Figure 4A force-displacement curve is shown. On curve marked values of equivalent force (green markers in Figure 4A) and values of applied force in the same displacements (red markers in Figure 4A). The relative error was from 3 to 15%. The region a crack appears in experiments matches the region in simulation. Of course, there were some deviations, which can be explained by not accurate boundary conditions. Despite the simple function for mechanical properties (binarization by threshold), adequate results were achieved. It can be explained by the stiffness matrix integration way. It proposes a method to consider the bone tissue structure. In this case, nodal displacement in element conditioned not only by mechanical properties but by material distribution. Previously it was shown that bone volume fraction and material distribution influence much greater on mechanical properties than other morphological variables (Maquer et al., 2015). To estimate the distribution approximation by tensor is popular. But in practice, it means that at the first step such tensor should be calculated, and then the tensors should be taken into account in mechanical model and simulation. The described method allows taking into account structural properties indirectly by stiffness matrix integration way. Add it costs less laborious.
More than that, the proposed method allows using part of the CT data space. In this case, the corresponding dimensions of the initial mesh should be set. From a practical point of view after the region of interest is chosen boundary conditions should be applied. So bone can be checked for strength. Typical loading cases can be simulated. The quality of the solution can be estimated by normalized error. And analyses can be automated in terms of the normalized error and the stress invariants, e.g., Von Mises stress.

CONCLUSION
The CT-based FEM method is introduced. A three-dimensional isoparametric linear element of a continuous medium, built on the basis of computed tomography data, presented. The general algorithm to get the FE model is presented. For this purpose computational and post-processing grids were introduced. To estimate the quality of simulation normalized error was used. Three-point bending was simulated by the method and results were compared with provided experiments. Received results illustrate the effectiveness of the method.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Care Committee of Kazan State Medical University.

AUTHOR CONTRIBUTIONS
OG, NK, PB, TB, and OS design method, realized it on programming language and wrote the manuscript. AF, MB, OS, FF, RI, ES, and TB conducted experimental data for numerical calculation and contributed to manuscript writing. All authors have read and agreed to the published version of the manuscript.

FUNDING
This work was supported by the Russian Science Foundation (RSF grant No. 18-75-10027).