In this study, we propose a deep learning model to predict the proximal femoral strength calculated from QCT-based FEA by integrating WGS features and DXA-derived imaging features. The studied cohort was acquired from the LOS (35, 36). The LOS cohort is an ongoing research dataset (>17,000 subjects accumulated so far) with recruitment starting in 2011, aimed at investigating both environmental and genetic risk factors for osteoporosis and other musculoskeletal diseases (37, 38). All participants signed an informed-consent document before any data collection, and the study was approved by the Tulane University Institutional Review Board.

Peak BMD achieved and remained relatively stable at ages 20-50 years is most powerful in predicting BMD and risk to osteoporotic fractures later in life due to the relatively stable physiological and hormone status during this age period (39–42). A 10% increase in peak BMD would delay the onset of osteoporosis by 13 years (40). In comparison, a 10% increase in the age of menopause, or a 10% reduction in age-related bone loss would only delay the onset of osteoporosis by 2 years (40).

Therefore, in this study, we focus on a cohort of 931 male subjects, aged 20-50, consisting of 345 African Americans and 586 Caucasians, with available QCT images, WGS and DXA-derived features. The basic demographic information for the enrolled subjects is shown in Table 1 .

The QCT scans (GE Discovery CT750 HD system; 2.5 mm-thick slices; pixel size, 0.695–0.986 mm; 512 × 512 matrix) were acquired at Tulane University Department of Radiology. For each QCT slice, contours of the left proximal femur were labeled by well-trained operators, in consultation with our experienced researcher (J.H.K.). We developed in-house software for automated annotation, containing a previously developed deep learning-based segmentation model (43) and a thresholding algorithm with edge tracing (44), in combination with manual visualization and correction. The deep learning model proposed by Zhao et al. (43) achieved a Dice similarity coefficient of 0.9888, indicating that only minor manual modifications are required for annotating new QCT images. Using the annotated contours, we then used linear and nonlinear FE models to estimate the strengths of the proximal femoral under two loading conditions, single-limb stance and loading from a fall onto the posterolateral aspect of the greater trochanter (19).

The nonlinear FEA models simulated mechanical testing of the femur in which displacement is incrementally applied to the femoral head (7, 16–19). The computed reaction force on the femoral head initially increases, reaches a peak value (the load capacity or fracture load), and then decreases. To achieve this mechanical behavior, the FEA models employ heterogeneous isotropic elastic moduli, yield strengths, and nonlinear post-yield properties. These properties are computed from the calibrated QCT density (ρ_CHA , g/cm³) of each voxel in an element, which are then used to compute the ash density (ρ_ASH , g/cm³) (ρ_ASH =0.0633 + 0.887 ρ_CHA ), and ρ_ASH is used to compute mechanical properties. Each linear hexahedral finite element measures 2.5 mm on a side and the mechanical properties of the element are computed by averaging the values of each property over all voxels in the element, while accounting for the volume fraction of each voxel within the element. Together, these mechanical properties describe an idealized density-dependent nonlinear stress-strain curve for each element (16–19). Material yield is defined to occur when the von Mises stress exceeds the yield strength of the element. After yield, the plastic flow was modeled assuming a plastic strain-rate vector normal to the von Mises yield surface and isotropic hardening/softening. Displacement is applied incrementally to the femoral head, and the reaction force on the femoral head is computed at each increment as the distal end of the model is fully constrained. For the fall models, the surface of the greater trochanter opposite the loaded surface of the femoral head was constrained in the direction of the displacements while allowing motion transversely. The nonlinear FEA-computed proximal femoral fracture load was defined as the maximum FEA-computed force on the femoral head, i.e., the load capacity.

For phenotypes, we calculated three proximal femoral strengths under the two loading conditions: LF, NLF and NLS. The LF represents the load at the onset of fracture (7). To determine the LF, the factor of safety (FOS) at the centroid of each finite element in the model is calculated as the ratio of the yield strength of the finite element to the von Mises stress at the centroid of the element. The LF is defined as the force applied to the femoral head when the FOS values of 15 contiguous non-surface elements are equal to or less than 1.0 (7, 19). The NLF and NLS were calculated using the above-described method and represent the load capacity of the proximal femur (the maximum force of the femoral head can support) (19). The basic statistical information for the calculated proximal femoral strengths is shown in Table 2 .

The WGS of the human peripheral blood DNA were performed with an average read depth of 22× using a BGISEQ-500 sequencer (BGI Americas Corporation, Cambridge, MA, USA) of 350 bp paired-end reads (38). The aligned and cleaned WGS data were mapped to the human reference genome (GRCh38/hg38) using Burrows-Wheeler Aligner software (45). This process followed the recommended best practices for variant analysis with the Genome Analysis Toolkit (GATK) to guarantee precise variant identification (46). The HaplotypeCaller tool within GATK was employed to identify genomic variations, and we further enhanced the reliability of our variant calls through the application of the variant quality score recalibration method (46).

There were a total of 10,623,292 single nucleotide polymorphisms (SNPs) in the cohort with 935 subjects. For quality control, we removed genetic variants with missing rates larger than 5%, Hardy-Weinberg equilibrium exact test p-values less than 10^-4, and minor allele frequency (MAF) less than 5%. Individuals with a missing rate larger than 20% were also excluded. Since subjects from two races were enrolled, principal component analysis (PCA) was applied to the genotypes and generated principal component scores (PCs) to perform population stratification or admixture (47). In addition, the age, weight, height, and first 10 PCs were used as covariates in GWAS (18, 35, 36).

The genome-wide association analyses were performed to test the association between each of three phenotypes and SNPs from WGS. Suppose that there are N subjects in the analyses. Let y_i be the value of the i-th subject for a phenotype and g_i be the genotype for the i-th subject, where g_i is the number of minor alleles that the subject carries at a SNP. We assume that there is a total of C covariates and the covariates for i-th subject are { v 1 ( i ) , v 2 ( i ) , ⋯ , v C ( i )   } . For each SNP, the linear regression model is used to examine the effect of a SNP on a phenotype, as defined in Eq. 1.

where β is the effect size of the SNP on the phenotype after adjusting for the covariates. The aim is to test the null hypothesis that the SNP is not associated with the phenotype, which is equivalent to test H ₀: β = 0. The score test statistic under this model is defined in Eq. 2.

where U = ∑ i = 1 N y ˜ i g ˜ i and V = 1 N ∑ i = 1 N y ˜ i 2 ∑ i = 1 N g ˜ i 2 , y ˜ i and g ˜ i are the adjusted phenotype and genotype for i-th subject for the covariates, indicating that y ˜ i is the residual of y_i under the linear regression model defined in Eq. 3.

where g ˜ i is the residual of g_i under the linear regression model defined in Eq. 4.

Under H ₀, T_score follows a standard normal distribution (48).

For each subject, aBMD (g/cm²) at various skeletal sites (lumbar spine, hip, forearm, and total body) and body composition (fat/lean mass) were measured using a Hologic Discovery-A DXA (Hologic Inc., USA) by trained and certified research staff at Tulane Center for Biomedical Informatics and Genomics. To ensure quality assurance, the machine was calibrated daily using a phantom scan. The accuracy of BMD measurement was assessed by the coefficient of variation for repeated measurements, which was approximately 1.9% for femoral neck BMD (38). In addition, all the DXA images have been reanalyzed using the TBS iNsight software (Medimaps Group, Geneva, Switzerland) to obtain trabecular bone score (TBS). As a result, 196 DXA-derived imaging features were obtained and used as the imaging features in this study.

For quality control purposes: the DXA machine was calibrated daily, and long-term precision was monitored by phantoms with a coefficient of variation ≤0.7% for spine aBMD and a coefficient of variation ≤1.0% for hip aBMD (49). Mechanical malfunction, radiation quality, absorption coefficient, and tissue-equivalent materials were also checked and calibrated before the aBMD examination on a daily basis. The radiologist was licensed in the State of Louisiana and registered through the American Registry of Radiologic Technologists. The detailed DXA-derived imaging features are shown in Table S1 .

For each Variational autoencoder (VAE), proposed by Kingma et al. (50), is a latent variable generative model which learns the deep representation of the input data. The goal of VAE is to maximize the marginal likelihood of the data (a.k.a evidence), which can be decomposed into a sum over marginal log-likelihoods of individual features, as illustrated in Eq. 5.

where x ^{( i )} is the feature vector for i-th subject in the dataset { x ( i ) } i = 1 N , N is the number of subjects, z is a random variable in the latent space, q_ϕ is the posterior approximation of z with the learnable parameters ϕ, p_θ is the ground truth posterior distribution of z with the intractable parameters θ, and D K L ( · ∥ · ) represents the Kullback–Leibler (KL) divergence between the approximated posterior distribution and the ground truth posterior distribution. Because of the non-negativity of the KL divergence, the log-likelihood log p θ ( x ( i ) ) ≥ ℒ ( θ , ϕ ; x ( i ) ) . If the approximated posterior distribution q ϕ ( z | x ( i ) ) is identical to the ground truth posterior distribution p θ ( z | x ( i ) ) , then the log p θ ( x ( i ) ) = ℒ ( θ , ϕ ; x ( i ) ) . Therefore, ℒ ( θ , ϕ ; x ( i ) ) is called the evidence lower bound (ELOB), which is defined by Eq. 6.

Thus, minimizing the KL divergence is equivalent to maximizing the ELOB. To train the model explicitly and implement the loss function in a closed form, we parameterize the q_ϕ as a multivariate normal distribution (multivariate Gaussian distribution) with an approximately diagonal variance-covariance matrix. Then the analytical solution for the KL divergence is shown in Eq. 7.

where D is the number of the latent variables extracted by the VAE, and μ d ( i ) and ( σ d ( i ) ) 2 are the approximate mean and variance of the posterior distribution of d-th latent variable for i-th subject.

We extend the VAE from single-view input into multi-view input fashion. Notably, as the fact that the product of Gaussian distributions is also a Gaussian distribution, we apply the PoE to generate the common latent space for the variation inference with an analytical solution. Suppose that under the multi-view setting, we have the data in M views, x ₁, x ₂, … x_M . For the data in m-th view (m=1,…, M), a nonlinear function implemented by a neural network is employed as the encoder, denoted as q ϕ m ( z m | x m ( i ) ) , where ϕ_m represents the learnable parameters of the nonlinear function for m-th view. For each encoder, we estimate the mean vector and the variance-covariance matrix of multivariate Gaussian distribution for the approximate posterior distribution, denoted as μ m ( i ) and Σ m ( i ) for i-th subject, and we assume μ m ( i ) ∈ ℝ D is a vector and Σ m ( i ) ∈ ℝ D × D is a diagonal matrix where D is the dimension of the latent space. In our implementation, we employ multi-layer perceptron (MLP) as the encoder. To guarantee the positivity of the covariance, the output of the MLP is denoted as the log Σ m ( i ) first and then is converted to Σ m ( i ) using the exponential function. Formally, the encoder is defined in Eq. 8.

where z_m is the latent variable extracted by m-th view with the dimension of D × 1. M L P m μ and M L P m Σ are the neural networks for calculating mean and covariance, respectively. Let T m ( i ) = ( Σ m ( i ) ) − 1 , then the multivariate Gaussian distribution for m-th view is rewritten as Eq. 9.

where Δ m ( i ) = − 1 2 ( μ m ( i ) ) T T m ( i ) μ m ( i ) − D 2 log 2 π + 1 2 log | T m ( i ) | .

A PoE models the target posterior distribution of the common latent variable from multi-view as the product of the individual posterior distribution of the latent variable from single-view. According to Eq. 9, Δ m ( i ) is not related to the latent variable z_m . Therefore, for the following analysis, Δ m ( i ) is considered as a constant. In our MVAE-PoE, the PoE generates the common latent variable z using Eq. 10.

That is, the multivariate Gaussian distribution of the common latent variable is defined by the product of the multivariate Gaussian distribution of the latent variable extracted by m-th view. According to (51), the approximated posterior distribution of the common latent variable, z, is shown in Eq. 11.

where μ z ( i ) and Σ z ( i ) are the mean vector and variance-covariance matrix of the approximated posterior distribution of common latent variable for i-th subject. To make the neural network differentiable, we adopt the reparameterization trick (50, 52) to reparametrize the mean vector and the diagonal variance-covariance matrix of the multivariate Gaussian distribution, as shown in Eq. 12.

where ϵ z ∼ N ( 0 , I ) and ⨀ indicates the element-wise product. For each view, we employed MLP layers as the decoder to restore the features, which is denoted as M L P m d e c for m-th view. The graphical architecture of the proposed MVAE-PoE model is shown in Figure 1 .

Using the multivariate Gaussian distribution, the ELOB for MVAE-PoE is derived in an explicit form, shown in Eq. 13.

The first term in the RHS of Eq. 13 is defined as the cross-entropy between the reconstructed data and the original input, and the second term in the RHS of Eq. 13 is the KL-divergence between the approximated posterior distribution and the true posterior distribution. The analytical form of the KL-divergence is the same as Eq. 3 since we employ the multivariate Gaussian distribution with an approximately diagonal variance-covariance as the ground truth. Thus, the close-form solution for the loss function is shown in Eq. 14.

where μ d ( i ) and ( σ d ( i ) ) 2 are the approximate mean and variance of the posterior distribution of d-th latent variable for i-th subject, and x ^ m ( i ) represents the reconstructed feature vector for i-th subject from m-th view.

20% of the subjects are randomly chosen as the test set, and the rest of the data are used as the training set. Predicting the proximal femoral strength is treated as a regression task. As shown in Figure 1 , a linear regression model is employed to predict the proximal femoral strengths using the extracted latent variables, z.

For model evaluation, mean absolute error (MAE), mean absolute percentage error (MAPE), root mean squared error (RMSE) and R ²-score are employed. The definitions of MAE, MAPE, RMSE and R ²-score are shown in Eqs. 11-14.

where y_i is the ground truth of the proximal femoral strength and y ^ i is the model prediction. A lower MAE/MAPE/RMSE and a higher R ² -score indicate better performance. According to Eqs. 15-18, 0 of MAE/MAPE/RMSE indicates the perfect match. According to Eq. 18, R ² -score ranges from -∞ to 1, where 1 indicates the perfect match.

Similar to (53), a leave-one-out technique is adopted to identify the feature significance in each view. A feature is significant if the performance of predicting proximal femoral strength decreases significantly when this feature is replaced by zero. By ranking the performance drops, the significance of the feature is obtained.

We performed GWAS analysis for testing the association between each of the three types of proximal femoral strengths, including linear fall fracture load (LF), nonlinear fall fracture load (NLF), and nonlinear stance fracture load (NLS), and each of the single nucleotide polymorphisms (SNPs) after quality control. Manhattan plots of these three types of proximal femoral strengths are depicted in Figure 2 . Since the sample size of the Louisiana Osteoporosis Study (LOS) cohort was relatively small for genetic association studies, we expanded our search space to look at a much wider landscape of associations by selecting top 256 SNPs with the lowest p-values to extract the WGS features that are associated with each phenotype. These identified SNPs were used as WGS features for the downstream task. Meanwhile, 196 DXA-derived imaging features were employed.

We trained and tested the MVAE-PoE model using our workstation with a NVIDIA RTX 3090 GPU and an Intel core I9 CPU. The designed models were implemented using TensorFlow 2.5. We performed the grid search to optimize hyperparameters, and the searching space included the number of MLP layers in both encoder and decoder: 1, 2, 3; the dimension of common latent space: 32, 48, 64, 128 or 256; and the number of hidden units for each MLP layer: 32, 48, 64, 128, or 256. Table 3 shows the best performance achieved using our proposed MVAE-PoE model for predicting three proximal femoral strengths. Also, we performed experiments using the different combinations of these three views to test the effectiveness of information fusion. We plotted the model prediction and the ground truth of the three FEA-computed proximal femoral strengths with the best performance in Figure 3 . In Figure 3 , each subject is represented by a blue dot. The vertical axis is the predicted strength. The red dashed line indicates a perfect match, and the green dashed line is the linear regression result of the prediction.

Integrating information from two views significantly improved the performance of predicting proximal femoral strength. According to Table 3 , the proposed MVAE-PoE model achieved its best performance for LF, NLF and NLS prediction using WGS features and DXA-derived imaging features. For example, the proposed model improved the R ²-score to 0.5569 compared with 0.4866 using DXA features alone for LF prediction.

DXA-derived imaging features are significantly more important than the WGS features in terms of the prediction performance. For LF, NLF and NLS, using DXA features alone, the designed models achieved the MAPEs of 18.04%, 6.84%, and 7.95%, respectively. Integrating DXA features with WGS features, the MAPEs were lowered by 1.62%, 0.55% and 0.46%, respectively. This finding was consistent with clinical practice that DXA-derived imaging features correlate with bone weakness and fragility fracture (10), with site-specific DXA explaining approximately 55% of the variability in predicting proximal femoral strengths (54).

Proximal femoral strengths depend on WGS features; however, using WGS features alone, the model does not generate satisfactory prediction results. For predicting LF, using only WGS features increased the MAPE from 18.04% to 24.47%. For predicting NLF, using only WGS increased the RMSE from 363.58 to 466.02.

We compared our MVAE-PoE model to other multi-view integration methods for prediction tasks. The tested models include:

For the above algorithms, a linear regression estimator was applied to perform the prediction task using the extracted latent variables. The overall performance comparison for predicting the three proximal femoral strengths is shown in Table 4 . For the compared algorithms, the grid search was also performed to find the best hyperparameters.

Compared to other multi-view information extraction models, the proposed MVAE-PoE achieved the best performance for predicting all types of proximal femoral strengths. The MCCA, KCCA, KGCCA and SCCA are four machine learning-based methods and the AAE is a deep learning-based method. For the MCCA, KCCA, KGCCA and SCCA, we trained these models with different dimensions of latent variables; for KCCA and KGCCA, we further tested the linear, polynomial and radial basis function (RBF) kernels. Even with tremendous hyperparameter fine-tuning, these machine learning-based methods didn’t generate better performance than the designed MVAE-PoE models. For the AAE model, we employed the same grid search settings. However, the achieved MAPEs were 20.68%, 7.22% and 10.27% for LF, NLF and NLS prediction, which indicated inferior performance than MVAE-PoE. The CCA-based methods have been commonly used in data fusion or integration; however, CCA-based methods treat the modalities as linearly and multivariately correlated without considering the direction of the linear relationship (60). In this study, we demonstrate that MVAE-PoE enables more useful and generalizable representations by capturing the abstract relationship between the views for downstream tasks such as prediction tasks.

Our model, which has shown excellent performance in predicting proximal femoral strength by integrating information from multiple views, holds promise for other radiogenomics data analysis problems, such as cancer prediction. By combining radiological imaging features with genomic data, our model can uncover valuable insights into the development, progression, and treatment response of diseases. Leveraging the power of our model, we believe it can effectively analyze radiogenomics data to enhance prediction accuracy and contribute to advancements in personalized medicine.

We applied the leave-one-out method to determine the feature importance. The leave-one-out indicated that we replaced one specific feature by zero for each subject in the test set when evaluating this feature, and the replaced features are named as zero-filled features. We compared the MAE changes between using the raw features and the zero-filled features. If the MAE between the GT and the model prediction increased significantly, then the evaluated feature was a significant feature. For each model, we listed the top 15 most important features in Figure 4 .

For the prediction of LF, according to Figure 4A , 10 of the top 15 most important features were WGS features, and 5 feature was a the DXA feature. The trochanter BMD (TROCH_BMD) was the most important DXA feature. Trochanteric BMD is associated with trochanteric fracture in the elderly and is among the best predictors of femoral strength (61). For NLF prediction, according to Figure 4B , one of the top 15 most significant features were DXA-derived features and the remaining 14 were WGS features. For NLS prediction, according to Figure 4C , 10 of the top 15 most significant features were DXA-derived features. This was consistent with the previous findings that DXA features explained approximately 55% of the proximal femoral strength while the proximal femoral strength was also influenced by genetics (25).

For each important WGS feature, we mapped the SNPs into the corresponding genes. The correspondingly associated clinical traits that were reported in the GWAS Catalog between each mapped gene and clinical traits are shown in Table 5 . For LF, According to the meta-analysis using 339,224 subjects from 125 subjects including African Americans and Caucasians, TSPAN12 show a positive correlation with BMI (62). COX6C and CAPG showed a strong correlation with body-shape indices on subjects from UK Biobank datasets (63). For NLF, the detected most important genes, CAPG, also showed a strong correlation with body fat distribution (64). For NLS, ERBB4 was associated with obesity on subjects from UK Biobank (65), which contained 339,244 individuals.

Our proposed research not only addressed issues related to multi-view information fusion, but also leveraged the value of widely used DXA with information provided by genetic markers for predicting proximal femoral strength. Therefore, this study has the potential to significantly impact both research and clinical practice. In the AGES-Reykjavik data set, Fleps et al. demonstrated that using FEA-computed hip fracture load to predict hip fracture was better than using total femoral aBMD only (22). In addition, genetic markers are important for identifying subjects at risk of hip fracture through effects on proximal femoral strength (25, 66, 10; 21). Due to their biological nature, genetic factors may also control nano-level bone mechanical properties and may further facilitate the prediction of bone strength and the assessment of hip fracture risk.

The most significant scientific impact of this study is the development and validation of the first comprehensive and accurate model for patient-specific assessment of predicting proximal femoral strength using multi-view information fusion by deep learning. Our multi-view deep learning-based model incorporates WGS features and DXA-derived imaging features, which are directly or indirectly related to proximal femoral strength and hip fracture. Deep learning-based techniques can automatically extract features and build accurate prediction models. Further, using the leave-one-out technique, the designed models are highly interpretable, leading to the identification of specific factors predictive of proximal femoral strengths.

The most practical clinical impact is the development and validation of an interpretable prediction model for proximal femoral strength using WGS features and DXA-derived image features, rather than using QCT. It is difficult to implement QCT-based femoral strength and hip fracture risk assessment in clinical practice due to the high radiation dosage and limited availability of QCT-based FEA. Lochmüller et al. suggested that clinical assessment of femoral fracture risk should preferably rely on femoral DXA (54). Our results suggest that there is a strong potential for using a combination of DXA and genetic markers to develop practical models for hip fracture risk assessment in the future.

Our study primarily focused on a population aged between 20 and 51 years. While this age group provides valuable insights into proximal femoral strength prediction using DXA and WGS, it is important to acknowledge that exclusion of older subjects may constrain the generalizability of our findings to the elderly population. Hence, future research endeavors should consider incorporating a more diverse age range to enhance the applicability of our predictive model across a broader population spectrum.