Optical bone densitometry robust to variation of soft tissue using machine learning techniques: validation by Monte Carlo simulation

Abstract. Significance: To achieve early detection of osteoporosis, a simple bone densitometry method using optics was proposed. However, individual differences in soft tissue structure and optical properties can cause errors in quantitative bone densitometry. Therefore, developing optical bone densitometry that is robust to soft tissue variations is important for the early detection of osteoporosis. Aim: The purpose of this study was to develop an optical bone densitometer that is insensitive to soft tissue, using Monte Carlo simulation and machine learning techniques, and to verify its feasibility. Approach: We propose a method to measure spatially resolved diffuse light from three directions of the biological tissue model and used machine learning techniques to predict bone density from these data. The three directions are backward, forward, and lateral to the direction of ballistic light irradiation. The method was validated using Monte Carlo simulations using synthetic biological tissue models with 1211 different random structural and optical properties. Results: The results were computed after a 10-fold cross-validation. From the simulated optical data, the machine learning model predicted bone density with a coefficient of determination of 0.760. Conclusions: The optical bone densitometry method proposed in this study was found to be robust against individual differences in soft tissue.

properties of the soft tissue were randomly determined to represent a population with sufficient variance. In the bone tissue, the three-dimensional (3D) trabecular pattern was generated using Alan Turing's reaction-diffusion model 22 because it is difficult to assign the trabecular bone a single optical property representative of a BMD value. The reaction-diffusion model is a widely used mathematical theory for describing the pattern formation process in biology, which is observed in bone during the early stages of calcification and development. 23,24 Because trabecular bone has a 3D nonlinear structure, the simulation was performed using a voxel-based Monte Carlo (VMC) model, which is a 3D extension of the Monte Carlo model for steady-state light transport in multilayered tissues (MCML) by Wang et al. 25,26 The purpose of this study is to develop and validate an ML prediction model with diffuse light as training data using a Monte Carlo simulation. The Monte Carlo model simulated spatially resolved steady-state diffuse light data acquired in three directions of the biological tissue model assuming an ultradistal radius. If this method is sufficiently accurate, it can be used for simple and safe bone densitometry, thus contributing to the early detection of osteoporosis.
2 Materials and Methods

Generation of Biological Tissue Model
The procedure for generating the biological tissue model is shown in Fig. 1. The biological tissue model consisted of the bone tissue, subcutaneous tissue, and dermis. The bone tissue was composed of the cortical bone, trabecular bone, and bone marrow. In this section, we describe the four steps for generating the biological tissue model.
In the first step, the structural properties of the biological tissue model were determined randomly. The structural properties of the biological tissue models are listed in Table 1. The thicknesses of the dermis and the subcutaneous tissue were determined using random numbers of uniform probability with a range of 1 to 2 mm 27 and 1 to 6 mm, 28 respectively. The structural   29 The values showing the range (number-number) have a uniform distribution of probability, and the mean ± standard deviation has a normal distribution. The three states of cortical bone thickness and BV/TV correspond to their respective values.
properties of the bone tissue were based on measurements of the ultradistal radius of women by Boutroy et al. 29 The cortical bone thickness (C.Th) and trabecular bone volume (BV) ratio [BV/tissue volume (TV)] vary at different stages of osteoporosis (osteoporosis, osteopenia, and healthy subjects). Therefore, the C.Th and BV/TV were determined using random numbers, assuming a normal distribution with different mean and variance values for each stage of osteoporosis. The distributions of C.Th and BV/TV were correlated with a Pearson correlation coefficient r of 0.54. The BMD of bone matrix (mBMD) assumed 1.2 g∕cm 3 of fully calcified bone.
In the second step, the trabecular pattern was generated using Alan Turing's 22 reactiondiffusion model. The modeling was based on a report by Miura and Maini. 30 The governing equation of the diffusion-reaction model is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 6 0 3 (1) This equation is called the activator-inhibitor system and is represented by the nonlinear reaction functions fðu; vÞ and gðu; vÞ, as well as diffusion terms, where u is the concentration of the activator and v is the inhibitor. The nonlinear terms f and g are E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 5 1 6 The diffusion term Δ is the Laplace operator, and du and dv are the diffusion coefficients. In this study, du and dv were set to 0.0002 and 0.01, respectively. The governing equation [Eq. (1)] was solved by an implicit scheme 31 based on the Fourier transform on a 0.025 grid in a 720 × 720 × 720 3D grid space, assuming periodic boundary conditions. The calculation was repeated 100 times with dt ¼ 1, and the solution was confirmed to converge. The initial states of u and v were assumed to be uniformly distributed random numbers within a range of AE0.5.
In the first step, the border between the trabecular bone and bone marrow space was defined by the threshold u th . First, a u th was applied to ten trabecular patterns, which were determined by randomly generating u with 10 different initial conditions. Accordingly, 10 trabecular bones were created with a similar BV/TV, in which u ≥ u th was defined as the trabecular bone and u < u th as marrow space. By repeating this process with 13 u th s, as shown in Fig. 2, the relationship between u th and BV/TV was obtained and is expressed as follows: Fig. 2 Relationship between BV/TV and u th , where u th is the threshold of the activation factor u. u th ¼ −7.662ðBV∕TVÞ 3 þ 1.645ðBV∕TVÞ 2 − 0.452BV∕TV þ 0.604: Equation (3) was fitted to a polynomial function using the least-squares method. Next, the voxel size u was defined. Because u is composed of voxels that do not have a size, we set the voxel size to 24.5 μm. This voxel size was similar to the resolution of the μCT image. The generated trabecular bone appeared to replicate the trabecular pattern observed in the real bone, as shown in Fig. 3. In addition, the generated trabecular bone was comparable to the ultradistal radius with respect to the trabecular number (Tb.N, mm −1 ), 29 fractal dimension, 32,33 and structure model index (SMI), 34,35 as shown in Table 2. These values quantitatively evaluate the trabecular bone shape. Because of the anisotropy of the trabecular pattern, the trabecular space (Tb.Sp) and variance (Tb.Sp SD), as well as the trabecular thickness (Tb.Th), were smaller than those of the ultradistal radius. In this study, this model was selected as the trabecular bone model because it can represent BMD changes in BV/TV, and the trabecular pattern resembles the shape of the real bone. Bone histomorphometry measurements were derived using TRI/3D-BON-FCS (Ratoc System Engineering Co. Ltd., Japan).
In the fourth step, the bone tissue, subcutaneous tissue, and dermis were defined, and a biological tissue model was generated, as shown in Fig. 4. The biological tissue model was assumed to be the ultradistal radius, and the bone axial direction was set as the y axis. First, the size of  bone tissue was determined. To ensure sufficient size in the bone axial direction, the generated trabecular bone was copied and joined in the y axis direction. From the cross-sectional area of the ultradistal radius, 29 a square with a size of 17.15 mm per side was defined as the outer bone surface (OBS) in the x-z plane (Fig. 4). Thus, the size of the bone tissue was 17.15 mm on the x-z axis and 35.28 mm on the y-axis. This size is comparable to the area of the trabecular bone in the distal radius. 29 Next, a cortical bone with a size of C.Th was defined from the OBS toward the center of the trabecular bone, and the bone tissue was generated. The aBMD of the bone tissue was defined using BV/TV, the size of the OBS in the x-axis direction (l obs , cm), and mBMD as follows: Finally, the subcutaneous tissue and dermis were generated in order from the OBS to the outward direction of the bone tissue, as shown in Fig. 4. All of the above processes were coded in Python 3.8.

Monte Carlo Simulation
Light transport in the synthetic biological tissue model was simulated using the Monte Carlo method. To deal with a tissue model containing trabecular bone with a 3D nonlinear pattern (Fig. 4), a VMC was built by extending MCML 25,26 to three dimensions. To represent individual differences, the optical properties of soft tissues were randomly determined. The VMC, as a voxel element, was defined as a hexahedron of size l v (24.5 μm). Voxels are assigned eigenvalues linked to optical properties and addresses indicating their location, and each voxel has a 3D space with the center coordinates as the origin. Photon packets were launched orthogonally onto the center of the biological tissue model, as shown in Fig. 4, which corresponds to a collimated infinitely narrow beam of photons. The weight of the initial photon was set to 1, and the specular reflection of light was treated in the same manner as in the MCML. Because the VMC has boundaries in six directions, the condition for a photon packet to hit with the voxel boundary is as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 1 3 5 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 9 2 s ¼ − where s is the propagation distance of the photon packet in one step 25 and d b is the distance to the voxel boundary. ξ is a uniformly distributed random number between 0 and 1, and μ a and μ s are the absorption and scattering coefficients, respectively. d b is calculated using the direction vector μ ðμ x ; μ y ; μ z Þ and position p ðp x ; p y ; p z Þ of a photon packet inside the voxel as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 6 8 7 where a photon packet hits the voxel boundary orthogonal to the axial direction, thus yielding a minimum value for d b (i.e., if the expression involving p x and μ x is the smallest compared with the other two, then the photon packet hits the boundary on the y-z plane). The positive or negative axial directions of the hitting photon packets can be determined from the direction vectors. When a photon packet hits a boundary, s is updated with E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 5 4 1 If there is a difference in the refractive index between the next and current voxels, whether the photon is internally reflected or transmitted is determined by the Fresnel equations and random numbers in the same manner as in MCML. For example, when a photon packet is transmitted, the photon position is updated from p i to p iþ1 and address a i ða x ; a y ; a z Þ is updated to a iþ1 ða xþ1 ; a y ; a z Þ for the transmission direction, as shown in Fig. 5. Then, when the photon packet hits the boundary again, the same operations are repeated until Eq. (5) becomes negative. If a photon packet escapes from the tissue model, its position, vector, and weight are preserved.
The optical properties of the biological tissue models are listed in Table 3. The range of optical properties was determined by referring to literature values, assuming measurements with ballistic light at a wavelength of 850 nm. μ a and μ s of the soft tissue were determined using the measurements of Simpson et al. 36 with uniform random numbers. The optical properties of the dermis measured by the authors included those of blacks and whites, and they assumed that the dermis and epidermis were combined. μ s of the bone was derived from the equation of Ugryumova et al. 16 Because mBMD is the wet density in Ugryumova's equation, μ s was calculated after converting mBMD to wet density from the data of Williams et al. 38 Thus, the conversion equation between mBMD and μ s is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 3 1 5 Most of the bone marrow in the radius is composed of adipose cells after the age of 20 to 30 years. 39 Therefore, the optical properties of the bone marrow are the same as those of the subcutaneous tissue, which is mainly composed of adipocytes.

Machine Learning
ML techniques were used to predict the aBMD from diffuse light simulated by the VMC. An overview of the system is presented in Fig. 6. In this section, we describe the four-step calculation procedure for predicting aBMD and the module that infers the function that relates simulated diffuse light to aBMD.

Prediction of bone density
In the first step, VMC was used to simulate light transport in a biological tissue model. Calculations were performed with a total photon count N of 10 7 . VMC was applied to 1211 tissue models with different structural and optical properties. In each model, the structural and optical properties were randomly combined in the ranges listed in Tables 1 and 3. Photon packets that escaped from the tissue model were categorized as backward, forward, or lateral to the direction of the packets launched. For the lateral direction, only the positive direction of the x-axis was considered. Photon packets escaping in three directions were defined as backscattered light (B), forward scattered light (F), and lateral scattered light (L).
In the second step, the physical quantity of packets escaping in the three directions was scored. For B and F, the sum of photon weights w was calculated for each radial distance r  from the photon packet launched coordinate with a range of Δr, as shown in Fig. 6. Here, r is represented by the array r ¼ ½0; Δr; 2Δr; 3Δr; : : : ; nΔr], and the sum of w is calculated from the range between r i and r iþ1 . The scored light intensity arrays B r and F r were E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 6 9 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 6 4 4 where I B and I F are the arrays that represent the sum of w per r i with respect to B and F, respectively, and ΔA r denotes the area of the annular ring, which is calculated as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 5 9 8 B r and F r were derived up to r ¼ 10 mm with Δr ¼ 0.4 mm. For B r , the initial position r 0 was set to 0.5 mm because the light intensity at r ¼ 0 was considered to be difficult to measure accurately when assuming actual measurements. For L, the sum of photon weights w was calculated for each positive distance along the z-axis from the photon packet irradiation coordinate with a range Δz, as shown in Fig. 6 (i.e., z is represented by the array z ¼ ½0; Δz; 2Δz; 3Δz; : : : ; nΔz, and the sum of w is calculated from the range bounded by z i and z i þ1 ). Here, the effective width of the y-axis is AEΔy. The scored light intensity array L z in the z-axis direction is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 4 6 9 where I L is an array that represents the sum of w per z i with respect to L and ΔA z indicates the area of the square, which is calculated as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 6 ; 4 0 2 ΔA z ¼ 2ΔyΔz: L z was derived up to z ¼ 20 mm with Δz ¼ 0.1 mm and Δy ¼ 5 mm.
In the third step, the obtained B r , F r , and L z were processed to generate the feature vectors. The feature vectors B f , F f , and L f are E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 6 ; 3 3 4 B f ¼ ½B r ; ln mB r ; ln vB r ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 2 9 0 F f ¼ ½ln mF r ; ln vF r ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 6 ; 2 6 7 L f ¼ ½ln L z ; ln mL z ; ln vL z ; where mB r , mF r , and mL r are the mean and vB r , vF r , and vL r are the variances of B r , F r , and L r , respectively. F r was not used in F f because F r did not form a valid distribution. For the L z signal, moving averages were calculated at 1 mm intervals. Then, to reduce the length of the L z array, averages of 2 mm intervals were adopted, that is, the length of array L z was 10. The elements of each feature vector were normalized to a mean of 0 and standard deviation (SD) of 1 for the dataset as follows: 40 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 1 1 6 ; 1 2 8 where μ and σ are the mean and SD of the dataset in each element of the feature vector and n is the total number of elements.
In the fourth step, aBMD was predicted using the ML model from B f , F f , and L f . The calculations with VMC, feature vector generation, and aBMD prediction using ML models were performed in Python 3.8.

Machine learning module
In this section, we describe the module that infers from labeled examples a function that relates the feature vectors (B f , F f , and L f ) and aBMD.
Four different ML techniques were tested: ridge regression (RR), support vector machine (SVM), random forest (RF), and gradient tree boosting (GTB). The criteria for selecting ML techniques were that they should be sufficiently stable and flexible for the particularities of the data, as well as have a strategy to control overgeneralization. RR solves a regression model in which the loss function is a linear least-squares function, and the regularization is given by the l2norm. 41,42 The coefficients of the regularization term were varied and tested in the range of 10 −5 to 10 −1 . SVM is less sensitive to noise using ε-insensitive loss functions and can construct nonlinear functions using kernel tricks. 43 A radial basis function (RBF) kernel was selected and tested with the kernel coefficient γ varying in the range of 10 −4 to 1. The soft margin C and tube ε were also adjusted. RF is one of the decision tree-based ensemble methods in which the output is an aggregation of the outputs of a set of classification and regression trees. 44 In RF, three parameters were adjusted: the number of decision trees, the minimum number of samples for a node to be considered a leaf, and the number of features to consider when calculating the optimal node split. 45 GTB is one of the decision-tree-based ensemble methods and is a generalization of boosting for arbitrary differentiable loss functions. 46,47 In GTB, six parameters were adjusted: the number of decision trees, the degeneracy of the step size used in the update to prevent overfitting, the maximum depth of the tree, the total minimum instance weights required for the children, the subsample ratio of the training instances, and the subsample ratio of the columns in building each tree. 45 For RR, SVM, and RF, we used the modules included in scikitlearn 0.24.2, and for GTB, we used XGboost 1.4.2.
To select the best structure and parameters for each ML module, 80% of the total data were used as a training and validation dataset (trDataset). The remaining 20% of the dataset (tsDataset) were used for testing and comparison purposes. The structure and parameter set of each algorithm were modified, and the best performing one was selected after a 10-fold cross-validation with grid search in trDataset. In 10-fold cross-validation, 90% of the dataset was randomly selected and used for training, and the remainder was used for validation. This was performed 10 times by rotating the dataset. The stability of the ML techniques was also verified by crossvalidation. After obtaining the best configuration, tsDataset was used to evaluate the performance. The coefficient of determination (r 2 ) was used as the metric for performance evaluation.

Results
The criterion for selecting the ML algorithm was the value of the coefficient of determination r 2 on tsDataset. tsDataset is a random selection of 20% (242) of the original dataset. The trDataset, which contains the remaining 80% (969), was used to train the algorithms. The ML algorithm was tuned in trDataset to achieve the best performance with 10-fold cross-validation. Once the best set of parameters and structures were found, the ML algorithms were retrained with the trDataset, and the performance was tested using tsDataset. Figure 7 shows a comparison of the coefficients of determination for each algorithm. The SVM regression exhibited the best performance (r 2 ¼ 0.757). The coefficients of determination for the other ML algorithms were r 2 ¼ 0.572 for RR, r 2 ¼ 0.451 for GTB, and r 2 ¼ 0.252 for RF. To determine the stability of the algorithm for unknown data, the SD of r 2 between the 10-fold cross-validation at trDataset was computed. The SD of r 2 in the cross-validation of SVM was 0.056, which is an acceptable value.
For the prediction of aBMD using SVM, the coefficient of determination r 2 on tsDataset was computed for all combinations of B f , L f , and F f (Fig. 8). The purpose of this study was to investigate the extent to which feature vectors and their combinations are related to aBMD. The parameters were tuned for all feature vector combinations with 10-fold cross-validation in trDataset. The prediction combining all of the feature vectors showed a coefficient of determination. Conversely, the prediction of aBMD from B f , L f , and F f alone exhibited lower performance (r 2 ≤ 0.302). These results suggest that it is difficult to predict aBMD with one feature vector, but combining feature vectors allows for highly accurate predictions.
As the final test, the performance of the system was tested on all 1211 data cases using the SVM method. SVM was selected because it gave the highest r 2 value for the aBMD prediction. The performance was assessed using 10-fold cross-validation on the entire dataset. The relationship between the predicted and reference values of aBMD is shown in Fig. 9. A linear regression of predicted and reference aBMD yielded an r 2 value of 0.760, indicating reasonable agreement. The Bland-Altman (BA) plot is shown in Fig. 10. The BA plot is a method used to check the agreement and systematic errors between two measurement methods. 48 The BA plot indicates a moderate correlation coefficient r of 0.22, with a slight proportional bias. This proportional bias  may not be a problem in practice. The mean difference between the predicted and reference values was 0.00, with no fixed bias. The limit of agreement for the predicted values was AE0.124 g∕cm 2 . These results suggest that BMD can be predicted with high accuracy using this method, even if there is variance in the thickness and optical properties of the soft tissue.

Discussion
The purpose of this study was to develop an optical bone densitometry method that is insensitive to individual differences in the structural and optical properties of soft tissues and to verify its feasibility. In the proposed method, spatially resolved diffuse light was obtained in three directions (backward, forward, and lateral to the direction of light irradiation) simulated by the Monte Carlo method, and BMD was predicted using ML techniques from these data. The crossvalidation demonstrated that the proposed method can predict BMD with high accuracy and low error in the simulation. The results suggest that the proposed method can be applied to optical bone densitometry, which is robust to variations in soft tissue.  Because the biological tissue model synthesized in this study is randomly constructed based on the range that a living organism can exhibit, the data obtained by the Monte Carlo simulation are considered to reflect a population with sufficient variance. The range of soft tissue properties was determined randomly based on measurements. In particular, the optical properties of the dermis cover a wide range of people, from colored individuals to Caucasians. 36 In bone tissue, a trabecular pattern was generated by the reaction-diffusion model because it is difficult to assign the trabecular bone a single scattering coefficient representative of a BMD value. The equation of Ugryumova et al. 16 yields negative scattering coefficients for a range of trabecular bone BMDs. In addition, Pifferi et al. found no age-related changes in the scattering coefficient in calcaneal measurements using TSR. 19 Overall, the unique and irregular shape of the trabecular pattern may lead to a different light scattering process compared with a homogeneous medium. The trabecular patterns generated in this study were similar to those of real bone in terms of appearance, BV/TV, and quantitative geometric structure. In addition, BMD was randomly adjusted from severe osteoporosis to the range of normal subjects. Moreover, the Monte Carlo method is the gold standard for modeling light transport in tissues. 49 Therefore, we consider that the diffuse light simulated in this study represents the structural and optical properties of the biological tissue for a population of sufficient variance.
The results shown in Fig. 7 seemingly contradict those reported by Chung et al. 14 Those authors demonstrated that the transmitted light intensity and aBMD are strongly correlated in the measurement of the ultradistal radius using an 850 nm wavelength light. However, in our results, the aBMD estimated only from transmitted light does not show a high coefficient of determination. This inconsistency seems to be due to the dispersion of the population. Chung et al. focused on a limited sample size of (10 participants). BMD measurements using only transmitted light are probably limited to populations with small variations in soft tissue composition. Nevertheless, F f , when combined with R, R f , and L f , clearly increased the coefficient of determination. This result suggests that the combination of light measured in different directions reduces the error from soft tissue with individual differences.
There is a possible nonlinear relationship between BMD and spatially resolved diffuse light observed in several directions. As shown in Fig. 8, SVM showed a higher prediction performance than RR. The RR is an algorithm based on linear multiple regression, 41,42 whereas SVM is an algorithm that can be applied to nonlinear data by mapping using nonlinear kernel functions, such as RBF. 43 In other words, the difference in the coefficient of determination between RR and SVM predictions of aBMD is probably due to the difference in the ability to support nonlinearity. GTB and RF are also nonlinear algorithms, but because they are decision-tree-based algorithms, they are probably not suitable for application to diffuse light, which shows continuous variation with distance from the light source. Therefore, when predicting BMD from our data, an algorithm for nonlinear data, such as SVM, is considered necessary.
The aBMD predicted using SVM from the Monte Carlo simulated data had a high coefficient of determination and lower error (Figs. 9 and 10). Thus, it has been demonstrated that our method can evaluate BMD with high prediction performance even when the bone is surrounded by soft tissue with individual differences in thickness and optical properties.
We acknowledge that there are several limitations to this study. First, the simulation was only a theoretical test. However, the simulations provide theoretical and clear insights into the different tissues that affect diffuse light. In addition, we believe that the combination of Monte Carlo simulation and ML techniques has several implications for the development of noninvasive medical measurements using the light diffusion theory that are beyond theoretical verification. An ML model built with sufficient variance and a large amount of data has excellent generalization performance; however, in the field of medical measurement, there are often ethical barriers and difficulties in data acquisition, such as a limited number of cases and invasive measurements. The simulation, which can generate almost unlimited data if computational resources are available, could offer a promising solution to such problems. Second, the VMC has boundaries in only six directions, which is an obvious limitation when compared with mesh-based methods 49-52 that can represent more free boundaries. However, inside a light-scattering medium, such as a biological tissue with a complex structure, light is averaged by scattering. Therefore, an overly accurate representation of the curvature of the trabecular structure is probably not practical. Third, a simple rectangular biological tissue model was adopted in this study. This model did not consider the heterogeneity in soft tissues, blood vessels, and complex bone structures; however, this information could be implemented in the model using CT and MRI techniques. However, the potential importance of this model is that it allows for a simple and targeted discussion of random changes in soft tissue thickness and optical properties in the validation of optical bone densitometry. This model might provide useful information about the interaction between the bone and soft tissue in measurements using diffuse light. All of the limitations mentioned here are attributed to the fact that it is unclear how much of the actual biological tissue should be assumed in the model to represent the light diffusion phenomenon and its output in vivo. It is also possible that simpler models represent substantial physical phenomena. Therefore, it is necessary to verify this by model experiments using phantoms and clinical trials.

Conclusion
In this study, we developed optical bone densitometry using Monte Carlo simulation and ML techniques and validated this method by cross-validation. From the results obtained, we concluded that the aBMD predicted from spatially resolved steady-state diffuse light data acquired in the backward, forward, and lateral planes of the biological tissue model was robust to differences in soft tissue layers thickness and optical properties.

Disclosures
The authors declare no financial or commercial conflicts of interest.