Analysis of Trabecular Bone Mechanics Using Machine Learning

Abstract

“Bone remodeling” is a dynamic process, and mutliphase analysis incorporated with the forecasting algorithm can help the biologists and orthopedics to interpret the laboratory generated results and to apply them in improving applications in the fields of “drug design, treatment, and therapy” of diseased bones. The metastasized bone microenvironment has always remained a challenging puzzle for the researchers. A multiphase computational model is interfaced with the artificial intelligence algorithm in a hybrid manner during this research. Trabecular surface remodeling is presented in this article, with the aid of video graphic footage, and the associated parametric thresholds are derived from artificial intelligence and clinical data.

Introduction

The bones in the skeletal system of vertebra are composed of two main classes: the cancellous bone and the cortical bone. Cortical bone is the hard tissue that is designed to resist well to tensile stress, on the outer surface of the femur. It overlies bone marrow and cancellous (trabecular) bone that occurs in the interior of the femur. The distribution of cortical and cancellous bone varies greatly between individual bones. Cancellous and cortical bones are structurally predominant in the neighborhood of joints and in the central sections of a femur away from the joints respectively. In addition, cortical and cancellous bone differ in their development, architecture, function, proximity to the bone marrow, blood supply, rapidity of turnover time, and magnitude of age-dependent changes, fractures, and classical immobilization-induced bone loss.

Clear understanding of the mechanochemistry of “bone osteoporosis” is highly desired,^1-4 whether it is a study devoted to explore a bone resorbing vicious cycle, where the cancerous cells play a complex role, ⁵ or it is a postmenopausal bone resorption. This mechanochemistry can be explored with the aid of the mechanosensory cells. The osteocytes are the principal mechanosensory cells of bone, these cells are believed to be activated by shear stress from fluid owing through the osteocyte canaliculi. It is believed that higher frequency and low amplitude strains can maintain bone as effectively as low frequency and high amplitude strains. ⁶ The mechanosensory behavior of osteocytes has remained a question of debate over the past two decades.^6-8 Through the experimental as well as the theoretical debates, the osteocytes have been declared as the sensory cells activated by the shear stress, since the osteocyte canaliculi is surrounded by fluid that exerts such shear stress. It is also believed that there is a network between osteocytes and the lining cells. Such network is electrically coupled and is a source of communication to control the bone’s structural adaptation as well as the homeostasis.

The continuum modeling of cortical and cancellous bone turnover initiated three decades ago.^9,10 In these studies, computational mechanics approach was used to test the influence of mechanical factors on bone remodeling, which could successfully predict changes in the external shape and apparent density of bone. Anisotropic continuum models have been proposed to describe the evolution of both apparent bone density and the orientation of trabecular architecture.^11,12 Even though the anisotropy is not explicitly considered in the rate equation, the osteocyte-regulated bone remodeling theory ¹³ could predict the emergence of the anisotropic trabecular microstructure.

The thinning and reduction of bone trabecula is one of the osteoporotic abnormalities. The animal experiments conducted by different research groups reported profound architectural changes in the loss of trabecular connectivity and the conversion of trabecular plates to rods.^14,15

The studies focusing on the osteoporosis supported the fact that whatever the cause is (ie, postmenopausal osteoporosis, vicious cycle of cancer or age), the trabeculae are reported to degenerate from plates to rod-like structure (although the quantification rate is still not clear).^16-18 In this article, we have considered the bone trabeculae finite element (FE) modeling and the detailed imaging footage of the morphological history of the trabecular plates and rods provided better understanding of the quantification rate. The trabecula is demonstrated in the previous studies to be a prominent contributor (in its plate-like form) in bone mechanics and the bone mechanical loading was mainly sustained by axially aligned trabeculae.^19,20 In the recent studies, a positive correlation was reported between axial bone volume fraction along imaging axes (aBV/TV) and corresponding axial elastic moduli along all directions. ²¹

Improved management is required for the patients requiring treatments that may diminish bone quality. This may be achieved by predicting the bone resorption rate and the important parametric values associated with it. The purpose of presenting a mathematical model in this study is to predict such parameters and to demonstrate the thinning of trabeculae with the aid of simulations.

Studies have strongly supported the hypothesis that mechanical strength depends on the bone mass as well as on the bone structure. ²² Therefore, clear understanding of the linkage between the bone trabeculae structure and the corresponding mechanical strength can address several questions and can help not only in better administration but also in proposing precautionary measures, for the patients.

In this article, we have proposed a novel hybrid approach; in a sense that, first machine learning is used to forecast the limiting values of the parameters; then these values are incorporated in the mathematical model to simulate the bone resorption. The architecture of the trabecula is directly modeled and the underlying objective of this article is (a) to understand the physical changes in trabecular structure associated with progression of bone resorption; (b) to evaluate the mechanical environment at the trabecular level, and (c) to propose a quantitative approach for trabecular surface remodeling. In the next section, the model is presented. The numerical results are presented in section “Results,” accompanied with the parametric values and graphical analysis; in the last section, some important conclusions are drawn.

Materials and Methods

Machine learning regression classifier

Over the past few decades, different variables governing the bone structural and mechanical properties are documented. These variables were extracted from the experimental and theoretical models. The experimental runs were concluded to be significant based on the traditional statistical methods; for example, certain statistical methods were used to obtain the correlation between the bone-age and Young’s modulus (“normalized stiffness”) as well as the failure-energy ²³ ; similarly, the regression analysis is used to explore the micro-CT images. ²⁴ With the advancement in research and technology, these statistical techniques have evolved over the past few years for better accuracy and for swift data handling. Different machine learning and artificial intelligence packages support the robust statistical tools, with the aid of which, the experimental data can be interpreted (by interfacing with the models), for better forecasting and clinical adminstration. ²⁵

The field of machine learning is aimed for the development of predictive models resulting in best possible accuracy and lowest error. Linear regression is considered for the understanding between input and output variables. ²⁶ However, regression classifier is aimed to discriminate the objects into classes. In a nutshell, it either predicts the true class of an object or classifies the data based on the training set.

In this article, the structural features of the bone trabecula are mimicked by the mathematical model and the morphological changes under the axisymmetric compressive loading are predicted using the machine learning approach. The clinical data ²⁷ for volume fraction as a response was interfaced into the model, while the yield stress and yield strain were taken as the predictor variables.

Model development

A trabecula model is developed with uniform distribution of $n$ sensor cells. Tube geometry is considered to better analyze the morphology of trabecula. The strain energy density per unit of mass can be described as

$$S_i=\frac{W_i}{{\rho }_i}$$

(1)

where $S_i$ , $W_i$ , and ${\rho }_i$ are the signal sensed by the sensor, the strain energy density, and the density at the ith sensor location. A single point $\mathbf{x}$ is expressed in its cylindrical coordinates ( $r$ , $\theta$ , $z$ ) and the geometry is illustrated in Figure 1. The density at location $\mathbf{x}$ is regulated by mean effect (stimulus) of the signals from all $m$ neighboring cells. The stimulus is mathematically expressed as

$$\psi (\mathbf{x},t)=\sum _{i=1}^{m}H(\mathbf{x})(S_i-S_{ref});0\le m\le n$$

(2)

and the substantial derivative of density as

$$\frac{D}{Dt}\rho (\mathbf{x},t)=C\psi (\mathbf{x},t)$$

(3)

where $H(\mathbf{x})=e^{-(\mathbf{x}-{\mathbf{x}}_{i2}/\mathbf{x}-{\mathbf{x}}_{k2})}0\le k\le n$ ; $C$ is the rate regulating time constant and ${\mathbf{x}}_k$ is defined as a location at which $H(x)\to e^{-1}$ . Equation (3) gives density at location $\mathbf{x}$ of the trabecula. This value must be less than the maximal density of trabecula ( $1.53\pm 0.25{\mathrm{g}/\mathrm{cm}}^3$ ).^28,29

Figure 1.

Schematic description of the two phases.

Naturally, the trabecula is surrounded by marrow; therefore, its rheological and mechanical properties must also be considered in the FE analysis; thus, the model of Mullender et al ¹³ is extended by taking into account the following system:

$$\frac{D}{Dt}{\rho }^{Mr}(\mathbf{x},t)=C^{Mr}{\psi }^{Mr}(\mathbf{x},t)$$

(4)

$$\frac{D}{Dt}{\rho }^{Tr}(\mathbf{x},t)=C^{Tr}{\psi }^{Tr}(\mathbf{x},t)$$

(5)

where the superscripts $Mr$ and $Tr$ stand for marrow and trabecula, respectively.

The mechanical stimulus for the marrow is taken into account in a manner similar to Birmingham et al, ²² whereas Equations (2) and (3) are considered for the trabecula.

To solve the system 4-5, we have considered a level set method, coupled with the axial symmetric stress-strain model for the stimulus and density evaluation. The stress-strain coupling with multiphase modeling has successfully been used to depict bone remodeling. ³⁰ The level set method has been used in the literature for analyzing the dynamics of elastic solids.^31-33 The advantage of this method is that the solver easily distinguishes the change in the rheological and mechanical properties of marrow and cancellous bone and thus helps depicting the resorption/formation caused by the mechanical stress.

The level set method is a smooth computational method, compared with front-tracking method,^34,35 which is one of the frequently used multiphase approaches of solid mechanics. The idea behind the level set method is to define a function $\varphi$ which is zero along the interface. The evolution equation for the interface is given by evolving $\varphi$ . When this function is re-initialized properly, it defines at all times a signed distance function to the interface. The numerical methodology for conservation laws captures the evolution of $\varphi$ .

The function $\varphi$ is defined as $1,-1$ , and $0$ for the trabecula, marrow, and interface, respectively. Figure 2 elaborates the phases in a schematic form. A Heaviside function defines the volume of fraction based on $\varphi$ values as

$$W(\varphi )=\{\begin{array}{ll}0\hfill & -1<\varphi <0\hfill \\ 1/2\hfill & \varphi =0\hfill \\ 1\hfill & 0<\varphi <1\hfill \end{array}$$

Figure 2.

Schematic of a single trabecula; the transverse view shows the surface cells (including osteoclasts and osteoblasts, canaliculi and osteocytes).

Therefore, the density equation is

$$\rho ={\rho }^{Mr}+({\rho }^{Tr}-{\rho }^{Mr})W(\varphi )$$

(6)

The trabecular architecture naturally adapts to the loading situation as per its location within the bone and on the external loading of the bone under consideration. In this model, the loading of the bone specimen has been mimicked from a general trabecula model, ²⁹ corresponding to a compression of the trabecular bone.

The coupled model and the corresponding conditions rely on the input values listed in Table 1. The stress, strain and the strain energy density are calculated in domains $A$ , $B$ , $C$ , $D,$ and $E$ (Figure 3) at each element using the FE analysis. The main assumption is that each element contains a sensor, and thus, the stimulus is calculated at each sensor location. The density is calculated and updated at each element using Equations (4) to (6). In case of no resorption/formation, the solver stops, if there is variation in trabecula structure and density, the solver updates itself and restarts (goes to the first step). The flow chart for the numerical solver is presented in Figure 4. We solved the system using a time-dependent solver. A parallel sparse direct and multi-recursive iterative linear solver has been used, which is robust and memory efficient. This method is used in the literature for solving large sparse symmetric, as well as nonsymmetric, linear equations.^37,38

Table 1.

Trabecula bone tissue and marrow material properties.

Poisson’s ratio $P^{Tr},P^{Mr}$	Young’s modulus $E^{Tr},E^{Mr}$
0.3, 0.167 - 36	$2915{({\rho }^{Tr})}^2,$ (0.008-0.02) MPa 7

Figure 3.

Longitudinal view of a single trabecula. $A$ and $B$ shows the marrow region, $C$ and $D$ represents the surface, $E$ represents the internal part, and $F$ shows the normal to the interface.

Figure 4.

Flow chart of the FE analysis. FE indicates finite element.

The simulations were extracted from the FE solver with relevant tolerance of 0.01, absolute tolerance of 0.001 and for a smaller time step, the time domain was scaled for simulations purpose and was then used for better understanding of the trabeculae resorption in days and years. The solver depicts resorption when mechanical properties of bone elements converge to those of marrow elements, and formation when the reverse process takes place. The dynamics of the interfacial elements undergoes an adaptive process of bone marrow or marrow bone and thus have material properties intermediary.

Results

In this article, the bone trabecula resorption is modeled. The parametric values are taken from an in vivo study. In Table 2, the parametric values are listed for the control and a group reporting bone resorption. We have used certain relationships ³⁹ to adapt the values in our numerical solver for a single trabecula. A step-by-step numerical solution was obtained. Two-dimensional stress analysis was analyzed by finite element modeling using diagonal scaling. ⁴⁰ Initially, the diameter of the trabecula was kept 12 mm. Von Mises equivalent stress was calculated for all the elements on trabecula surface using the relation for normal stress-strain for the coupled system. Trabecula/marrow interfacial surface movement was depicted using Equations (4) to (6) and Table 2. Convergence was achieved for 11 253 degrees of freedom, and for 976 number of mesh points. The minimum element quality was 0.1859 and element area ratio was 0.25.

Table 2.

Volumetric density and trabecular microarchitecture in a control group and a group with osteoporosis.

	Control group	Osteoporosis group
Area (mm2)	290.2 ± 67.2	291.7 ± 60.8
vBMD (mg HA a /cm3)	170.0 ± 27.0	163.3 ± 36.6
BV/TV (%)	0.154 ± 0.024	0.143 ± 0.026

a

Hydroxyapatite.

BV/TV, bone volume to trabacula volume ratio; vBMD, volumetric bone mas density.

In Figure 5, bone resorption is depicted (white region shows the deformed interface, gray region shows the marrow, and black shows the trabecula). The arrows are pointed toward the regions where the osteoblasts activity is influenced, resulting in resorption. In Figures 6 and 7, the resorption with the passage of time is presented as a result of the therapy, which causes an imbalance in the osteoclast-osteoblast activity, and these effects are adapted in the model through the specified boundary conditions along the interface.

Figure 5.

Bone resorption mechanism, white region shows the modified interface, gray region shows the marrow, and black shows the trabecula, the arrows shows direction the interfacial normal vector.

Figure 6.

With the passage of time, the osteoclastic resorption process grows abnormally deep and severs the trabecula.

Figure 7.

Wired network of cells in marrow (pink), trabecula (blue), and internal trabecula (dark blue).

Conclusions

During this study, we demonstrate that marrow/trabecula interfacial dynamics can be easily depicted using a computational model. The mechanical stress and the density gradients due to the signaling mechanism (between the osteocyte, osteoblast and osteoclast) can be well demonstrated using the coupled level set FE analysis. We have showed with the aid of our numerical results, that during osteoporosis, the original trabecular bone gradually transforms from an intact to more porous and damaged structure during the simulation of bone remodeling.

Our future work focuses on building anisotropic poro-elastic FE models that include the vascular porosity of the trabecula. This will help to provide a clear picture of the bone resorption dynamics during secondary osteoporosis, to compare with the experimental findings and to make future predictions.