WO2007025155A2 - Systems, products, and methods for predicting changes and fracture in trabecular bone - Google Patents

Abstract

Systems, products, and methods for predicting changes and fracture in trabecular bone are described herein. In some embodiments, systems, products, and methods for simulating bone loss that include reading a three-dimensional image of a trabecular bone, determining a number and location of bone resorption sites, and carrying out simulated bone resorption are described. In other embodiments, systems, products, and methods for determining risk of fracture in trabecular bone that include reading a three-dimensional image of a trabecular bone; removing voxels from trabeculae portion of the image to obtain a skeletonized structure of the trabecular bone, where the skeletonized structure of the trabecular bone is represented by voxels classified according to at least one topological classification; and adding reconstruction voxels to at least voxels of a first topological classification are described.

Description

SYSTEMS, PRODUCTS, AND METHODS FOR PREDICTING CHANGES AND FRACTURE IN TRABECULAR BONE

CROSS-REFERENCE TO RELATED APPLIATIONS

[0001] This application claims the benefit of the filing date of U.S. Patent Application No. 60/711,059, filed on August 24, 2005, and U.S. Patent Application No. 60/749,696, filed on December 12, 2005, the contents of which are hereby incorporated by reference herein in their entireties.

TECHNICAL FIELD

[0002] The disclosed subject matter relates to systems, products, and methods for predicting changes and fracture in trabecular bone. For example, the disclosed subject matter relates to systems, products, and methods for predicting morphological changes of microstructure and mechanical properties such as elastic modulus and strength in trabecular bone.

BACKGROUND

[0003] Generally speaking, sections of bone can be categorized into two groups: trabecular bone and cortical bone. Trabecular bone is the sponge-like bone in the ends of long bones and vertebrae and cortical bone is the hard outer shell of bones. Trabecular bone is typically more prone to a bone remodeling process that occurs naturally where bones are constantly broken down (also called resorption) and rebuilt (also called formation). In humans, osteoclasts break down the bone, converting the calcium salts in the bone to a form that passes easily into another form, while osteoblasts produce organic fibers on which calcium salts are deposited. In normal circumstances, the activities of these cells (osteoclasts and osteoblasts) are balanced as to maintain a constant bone mass in the body. [0004] Microstructure of trabecular bone consists of three-dimensional network of individual bone units called trabecula(e). These individual trabeculae have in general two types of morphology: rod-like and plate-like. The mechanical properties of trabecular bone depend on heavily on the amount of rod-like or plate-like trabeculae, the connections between rod-like and rod-like, or rod-like and plate-like, or plate-like and plate-like trabeculae, orientations of rod-like or plate-trabeculae, and the dimensions of rod-like or plate-like trabeculae. In osteoporosis, there is alterations in microstructure of trabecular bone such as transition from predominate plate-like trabeculae to rod-like trabeculae, loss of number of trabeculae or thinning of trabeculae.

[0005] The remodeling process provides two purposes. First, remodeling allows reshaping of the bones in response to the various different local stresses applied to the bone to ensure enough support is provided where needed. Second, remodeling helps the calcium levels in the bloodstream to be maintained. Calcium serves many important functions in the body, such as muscle contraction, nerve function, blood coagulation, and cell division, and if the calcium levels in the bloodstream drop, more calcium is released from the bone into the bloodstream. Hence, when the net release of calcium from the bone into the bloodstream is greater than the net calcification onto the bone, overall loss of bone mass can occur. [0006] Research during the past three decades has suggested a correlation between bone density and the elastic modulus or strength of trabecular bone. It is generally believed that fracture risk is enhanced with decreasing bone density as rate of resorption exceeds the rate of formation. For example, osteoporosis is a disease that causes perforations within the trabecular bone, which can grow over time, leading to decreased bone mass density. Women are generally more susceptible to osteoporosis, especially when estrogen levels fall during and/or after menopause. Due to the reduced level of estrogen during and/or after menopause, an increased activity in the bone remodeling process, and a subsequent decrease in bone mass, can occur despite bone formation activities. In addition, bone density itself may not be the only factor in determining bone failure. For example, even at a given bone density, risk of fracture risk is believed to increase with increasing age.

[0007] It should be noted that predictive models of bone based on in vitro or in vivo high resolution image can be quite difficult or computationally expensive to implement in the clinical setting.

[0008] Therefore, there is a need to predict the fracture risks in trabecular bones before the onset of osteoporosis has taken place as a function of trabecular bone microarchitecture. Moreover, there is a need for computationally less expensive model for clinical predictions of fracture risk. SUMMARY

[0009] The disclosed subject matter relates to methods, products, and systems for predicting change and failure of bones.

[0010] In some embodiments, systems, products, and methods for simulating bone loss that include reading a three-dimensional image of a trabecular bone; determining a number and location of bone resorption sites; and carrying out simulated bone resorption are described.

[0011] In other embodiments, systems, products, and methods for determining risk of fracture in trabecular bone that include reading a three-dimensional image of a trabecular bone; removing voxels from trabeculae portion of the image to obtain a skeletonized structure of the trabecular bone, where the skeletonized structure of the trabecular bone is represented by voxels classified according to at least one topological classification; and adding reconstruction voxels to at least voxels of a first topological classification are described. [0012] From this, a morphological analysis of trabecular bone microstructure based on individual trabecular rods and plates can be performed. In addition, a computational model based on individual trabecular rods and plates can be used for determining mechanical properties.

BRIEF DESCRIPTION OF THE DRAWINGS

[0013] The above and other objects and advantages of the disclosed subject matter will be apparent upon consideration of the following detailed description, taken in conjunction with the accompanying drawings, in which:

[0014] FIG. 1 is a diagram showing methods for simulating osteoporosis and risk of failure of trabecular bones in accordance with certain embodiments of the disclosed subject matter;

[0015] FIG. 2 is an image showing the conversion of a bone image into an eight-node brick element in accordance with certain embodiments of the disclosed subject matter;

[0016] FIG. 3A is an image showing bone loss due to perforation in accordance with certain embodiments of the disclosed subject matter;

[0017] FIG. 3B is an image showing bone loss due to broken trabeculae in accordance with certain embodiments of the disclosed subject matter;

[0018] FIG. 3C is an image showing bone loss due to isolated trabeculae in accordance with certain embodiments of the disclosed subject matter;

[0019] FIG. 4 is a diagram showing methods for carrying out morphological and finite element analyses of trabecular bone in accordance with certain embodiments of the disclosed subject matter;

[0020] FIG. 5A shows an image of a bone obtained using μCT in accordance with certain embodiments of the disclosed subject matter;

[0021] FIG. 5B shows a skeletonized image of the bone of FIG. 5 A in accordance with certain embodiments of the disclosed subject matter;

[0022] FIG. 6A shows an image of a bone having a perforation obtained using μCT in accordance with certain embodiments of the disclosed subject matter;

[0023] FIG. 6B shows a skeletonized image of the bone of FIG. 6 A maintaining the perforation in accordance with certain embodiments of the disclosed subject matter;

[0024] FIG. 7A shows an image of a bone obtained using μCT in accordance with certain embodiments of the disclosed subject matter;

[0025] FIG. 7B shows a skeletonized image of the bone of FIG. 7A in accordance with certain embodiments of the disclosed subject matter;

[0026] FIG. 7C shows an arc skeletonized image of the bone of FIG. 7 A in accordance with certain embodiments of the disclosed subject matter; [0027] FIG. 7D shows a reconstructed image of the bone of FIG. 7 A in accordance with certain embodiments of the disclosed subject matter;

[0028] FIG. 8A shows an image of a bone obtained using μCT in accordance with certain embodiments of the disclosed subject matter;

[0029] FIG. 8B shows a plate-reconstructed image of the bone of FIG. 8A in accordance with certain embodiments of the disclosed subject matter;

[0030] FIG. 8C shows a rod-reconstructed image of the bone of FIG. 8 A in accordance with certain embodiments of the disclosed subject matter;

[0031] FIG. 9 shows an image of calculating of average geometric vectors for trabecular plates and rods in accordance with certain embodiments of the disclosed subject matter;

[0032] FIG. 10 shows a schematic representation of trabecular meshing utilizing a

Delaunay triangulation technique in accordance with certain embodiments of the disclosed subject matter;

[0033] FIG. HA shows an image of a bone obtained using μCT in accordance with certain embodiments of the disclosed subject matter;

[0034] FIG. HB shows a skeletonized image of the bone of FIG. HA in accordance with certain embodiments of the disclosed subject matter;

[0035] FIG. HC shows a skeletonized image of the bone of FIG. HA with the voxels classified as a surface, curve, edge, and junction, in accordance with certain embodiments of the disclosed subject matter;

[0036] FIG. HD shows a fully reconstructed image of the bone of FIG. HA in accordance with certain embodiments of the disclosed subject matter;

[0037] FIG. HE shows a rod-reconstructed image of the bone of FIG. HA in accordance with certain embodiments of the disclosed subject matter;

[0038] FIG. HF shows a plate-reconstructed image of the bone of FIG. HA in accordance with certain embodiments of the disclosed subject matter;

[0039] FIG. 12 shows a graph showing the results of linear correlation analyses between the elastic modulus of the full voxel model and that of each of skeleton, plate-reconstructed, and rod-reconstructed models in accordance with certain embodiments of the disclosed subject matter;

[0040] FIG. 13 shows a graph showing the results of linear correlation analyses between the bone volume fraction of the full voxel models and that of each of skeleton, plate- reconstructed, and rod-reconstructed models in accordance with certain embodiments of the disclosed subject matter;

[0041] FIG. 14A shows images of full voxel and skeleton models of a trabecular bone having a plate-to-rod ratio of 20.10, plate tissue fraction of 0.9526, and a Young's modulus of 974.8 MPa in accordance with certain embodiments of the disclosed subject matter; [0042] FIG. 14B shows images of full voxel and skeleton models of a trabecular bone having a plate-to-rod ratio of 10.22, plate tissue fraction of 0.9109, and a Young's modulus of 674.8 MPa in accordance with certain embodiments of the disclosed subject matter; [0043] FIG. 14C shows images of full voxel and skeleton models of a trabecular bone having a plate-to-rod ratio of 5.94, plate tissue fraction of 0.8559, and a Young's modulus of 216.8 MPa in accordance with certain embodiments of the disclosed subject matter; [0044] FIG. 15A shows a graph showing the results of correlation analysis between Young's Modulus of the full voxel model and that of plate bone volume fraction by nonlinear regression of power laws in accordance with certain embodiments of the disclosed subject matter;

[0045] FIG. 15B shows a graph showing the results of correlation analysis between Young's Modulus of the full voxel model and that of rod bone volume fraction by nonlinear regression of power laws in accordance with certain embodiments of the disclosed subject matter;

[0046] FIG. 15C shows a graph showing the results of correlation analysis between Young's Modulus of the full voxel model and that of plate tissue fraction by nonlinear regression of power laws in accordance with certain embodiments of the disclosed subject matter;

[0047] FIG. 15D shows a graph showing the results of correlation analysis between Young's Modulus of the full voxel model and that of rod tissue fraction by nonlinear regression of power laws in accordance with certain embodiments of the disclosed subject matter;

[0048] FIG. 15E shows a graph showing the results of correlation analysis between Young's Modulus of the full voxel model and that of bone volume fraction by nonlinear regression of power laws in accordance with certain embodiments of the disclosed subject matter;

[0049] FIG. 15F shows a graph showing the results of correlation analysis between Young's Modulus of the full voxel model and that of plate-to-rod ratio by nonlinear regression of power laws in accordance with certain embodiments of the disclosed subject matter;

[0050] FIGS. 16A and 16F show images of bone obtained using μCT in accordance with certain embodiments of the disclosed subject matter;

[0051] FIGS. 16B and 16G show skeletonized images of the bone of FIGS. 16A and 16F, respectively, in accordance with certain embodiments of the disclosed subject matter;

[0052] FIGS. 16C and 16H show arc skeletonized images of the bone of FIGS. 16A and

16F, respectively, with the voxels classified as a surface, curve, edge, and junction in accordance with certain embodiments of the disclosed subject matter;

[0053] FIGS. 16D and 161 show partially reconstructed images of the bone of FIGS.

16A and 16F, respectively, in accordance with certain embodiments of the disclosed subject matter;

[0054] FIGS. 16E and 16J show fully reconstructed images of the bone of FIGS. 16A and 16F, respectively, in accordance with certain embodiments of the disclosed subject matter;

[0055] FIG. 17A shows images of isolated plate trabecula under different amounts of applied strain in accordance with certain embodiments of the disclosed subject matter;

[0056] FIG. 17B shows images of isolated rod trabecula under different amounts of applied strain in accordance with certain embodiments of the disclosed subject matter;

[0057] FIG. 18 shows a number of yielded trabecular plates and rods with increasing strain in accordance with certain embodiments of the disclosed subject matter;

[0058] FIG. 19A shows a graph of number of yielded trabecular plates in horizontal, oblique, and longitudinal orientations with increasing strain in accordance with certain embodiments of the disclosed subject matter;

[0059] FIG. 19B shows a graph of number of yielded trabecular rods in horizontal, oblique, and longitudinal orientations with increasing strain in accordance with certain embodiments of the disclosed subject matter;

[0060] FIGS. 2OA, 2OB, and 2OC show images of vertical, oblique, and horizontal rods, respectively, segmented from 3D μCT images of trabecular bone in accordance with certain embodiments of the disclosed subject matter;

[0061] FIG. 21 shows a graph showing the changes in the bone volume fraction, the

Young's modulus along the x (E_xx), y (E_yy), and z (E_zz) direction, and shear modulus along the xy (G_xy), yz (G_yz), and xz (G_xz) directions in accordance with certain embodiments of the disclosed subject matter;

[0062] FIGS. 22A, 22B, and 22C show bone loss due to perforation, broken trabeculae, and isolated trabeculae, respectively, in accordance with certain embodiments of the disclosed subject matter;

[0063] FIG. 23 shows a graph of the activation frequency change during menopause in accordance with certain embodiments of the disclosed subject matter;

[0064] FIGS. 24A, 24B, 24C, and 24D show simulations of a trabecular bone sample at 5 years before menopause, 0 years after menopause, 5 years after menopause, and 10 years after menopause, respectively, with perforations shown in the insets and breakage and disconnections shown in circles in accordance with certain embodiments of the disclosed subject matter;

[0065] FIG. 25 shows a graph showing the averaged, relative bone volume fraction due to various types of bone loss during menopause and the averaged, relative elastic modulus during menopause in accordance with certain embodiments of the disclosed subject matter;

[0066] FIGS. 26A and 26B show a graph showing the predicted femoral neck and spine bone volume fraction change compared with the corresponding clinical data in accordance with certain embodiments of the disclosed subject matter;

[0067] FIG. 27A shows a graph of the time course of averaged relative bone mass density from simulations of eight samples during treatment of bisphosphonate due to one or combination of mechanisms: reduced activation frequency (RS); increased secondary mineralization (MIN); and positive bone balance (PBB), compared with clinical data in accordance with certain embodiments of the disclosed subject matter;

[0068] FIGS. 27B and 27C show 2D distributions of trabecular bone with secondary mineralization before and after treatment, respectively, in accordance with certain embodiments of the disclosed subject matter; and

[0069] FIG. 28 shows a system that includes a processor, a memory, a display device, an input device, and an optional interface to an image capturing device in accordance with certain embodiments of the disclosed subject matter. DETAILED DESCRIPTION

[0070] As shown in FIG. 1, methods, products, and systems in accordance with certain embodiments of the disclosed subject matter can be used to predict change and fracture in trabecular bone. As shown in FIG. 1, this can include generating or obtaining an image of bone (100 in FIG. 1). The image can be a digital image. The image can be a three- dimensional (3D) image, hi certain embodiments, the image can be obtained by micro- computed tomography (μCT) or micro magnetic resonance images (μMRI). The obtained images can be represented as voxels (i.e., a 3D equivalent of a pixel) based on the sample size and the scanning resolution of the measurement. For example, a 4x4x4 mm cubical sample scanned at 21 μm nominal isotropic resolution can have approximately 191x191x191

Suitable thresholding, such as adaptive thresholding, global

thresholding, and bone volume fraction ("BVF") compensated global (or "matched global") thresholding, can be carried out. Such thresholding may be able to account for shading, noise, and the like that may have been introduced during the sample measurement process. In adaptive thresholding, local (adaptive) thresholds can be obtained in local volumes throughout the 3D sample depending on the greyscale value of the surrounding voxels. In global thresholding, a unique thresholding value can be selected based on the historgram distribution of the grey level values in the 3D image. In matched global thresholding, global threshold value can be adjusted so that the BVF approximately matches the BVF obtained through physical density measurements. Any other suitable image manipulation techniques can be employed as will be readily apparent to one of ordinary skill in the art. [0071] As an example of obtaining an image of a bone, vivaCT 40 (SCANCO Medical AG, Switzerland) can be utilized to capture μCT images using an isotropic, nominal resolution of 21 microns. The obtained image may be a 16-bit image and can be converted to an 8-bit image. The reconstructed μCT image can be filtered to remove noise, and the cleaned μCT images can then be thresholded to distinguish between bone voxels and bone marrow voxels using a global value thresholding technique. Then, isolated voxels or voxel- clusters can be removed using a clustering analysis for further analysis, such as morphological and finite element analyses (e.g., 102 and 104 of FIG. 1), as will be described below. Other mechanisms for capturing μCT images can also be used. [0072] For an example of μMRl image collection, bone marrow contained in the bone can be removed to prevent image artifacts that may arise from air leaking into the marrow spaces and from paramagnetic deoxyhemoglobin in erythrocytes of hematopoietic marrow. μMRI images can be acquired on a Bruker Instruments DMX-400 spectrometer/micro- imaging system operating at 400 MHz proton resonance frequency (9.4 Tesla magnetic field strength). This system can be fitted with a micro-imaging accessory, consisting of 3-axes magnetic field gradients, with 100 Gauss/cm maximum gradient amplitude in all three channels. A 1.5 cm diameter RF coil can be used to obtain high detection sensitivity with sufficient RF field homogeneity. For the pulse sequence, a partial flip-angle gradient-echo technique (FLASH-type sequences, Haase) can be utilized. This pulse sequence can be tailored to be insensitive to susceptibility-induced intravoxel phase dispersion at the bone- marrow interface (a phenomenon which is known to cause local signal loss and leading to an apparent thickening of trabeculae in FLASH-type sequences). Images can be acquired with an echo time (TE) of 8 msec, one signal average, and pulse repetition time (TR) of 60 msec, or any other suitable technique to provide an adequate signal to noise ratio (SNR). The optimum radio frequency (RF) flip angle α under these conditions can be 140 degrees. The field of view can be (10mm)³, corresponding to an isotropic linear resolution of 68.5x68.5x102.5 μm³. Other mechanisms for capturing μMRI images can also be used. [0073] Having obtained a suitable image of trabecular bone, finite element analysis (FEA) can optionally be performed to determine the stress field applied to the bone (106 of FIG. 1). For example, each bone voxel can be converted to an eight-node brick element image (see FIG. 2). Other suitable conversion of data to optimize computational time can be employed as will be readily apparent to one of ordinary skill in the art. Trabecular bone tissue can be modeled as an elastic-plastic strain-hardening material with an asymmetric tensile-compressive failure criterion, and a geometrically non-linear large deformation kinematic behavior can be assumed. Tension-compression asymmetry can be introduced by shifting the yield surface using pseudo kinematic hardening.

[0074] A parallel finite element software, OLYMPUS, can be utilized. Olympus, built on a serial finite element analysis program, uses a parallel graph partitioner, PARMETIS, and an unstructured parallel multigrid equation solver, PROMETHEUS. Constant material properties can be employed in the computation, although specimen specific calibrations for nonlinear analyses can also be utilized if desired. The predictions in the apparent elastic modulus, yield strength, and yield strain by voxel based FEA models can be compared with the experimental measurements.

[0075] All FEA models can be analyzed using custom software on an IBM-Power 4 (IBM Corp., White Plains, NY) distributed memory parallel supercomputer (Datastar located at San Diego Super Computing Center SDSC) available through the National Partnership for Advanced Computational Infrastructure (NPACI). Each analysis can be a large deformation and materially non-linear analysis. The apparent elastic modulus can be calculated from the initial linear step of the computation and the yield stress and strain can be determined by the 0.2% offset method.

[0076] As another example, FEA can be performed using a commercially available finite element software ABAQUS, as will be described below in greater detail. Any other suitable finite element analysis technique can be utilized, as will be readily apparent to one of ordinary skill in the art.

[0077] Subsequently, an activation frequency of a bone remodeling process, which can be obtained experimentally or from literature, can be input as a parameter (108 of FIG. 1). The activation frequency is a measure of how often a bone remodeling cycle (i.e., resorption and formation) occurs in the bone. It should be noted that the activation frequency need not be constant for each iteration, but can vary for each iteration. For example, the activation frequency can change very little before menopause, increase rapidly after the onset of menopause, and plateau off a few years after menopause. An another example, during treatment of osteoporosis, activation frequency can show a lower value, or even a negative value for successive iterations.

[0078] Based on the activation frequency parameter, areas where resorption can occur on the bone surface can be labeled and created as shown in 110 of FIG. 1. The resorption sites can be specified in any suitable manner. For example, the total number of resorption sites can be determined by the activation frequency and randomly distributed over the entire bone surface. As another example, voxels that have the highest stress fields, as determined in 106 of FIG. 1, can be labeled as resorption sites based on the activation frequency parameter. Other suitable methods for labeling the bone image with resorption sites can be employed, as will be readily apparent to one of ordinary skill in the art.

[0079] The resorption period, the resorption rate, and the wall thickness can also be provided as process parameters (112 of FIG. 1). The resorption period is the time it takes for the resorption process to occur before bone formation occurs, the resorption rate is the rate of bone loss near the resorption sites, and the wall thickness is the amount of bone lost in the resorption sites and can be related to the number and activity of osteoclasts at the resorption cavity. These process parameters can be obtained experimentally or obtained from literature. [0080] Based on the initial bone morphology, activation frequency, resorption period, resorption rate, and trabecula thickness described above, the resoprtion process can be carried out (114 of FIG. 1) near the resorption sites that were labeled in 110 of FIG. 1. As resorption proceeds, various different types of resorption cavities can be created. For example, resorption can lead to the formation of perforations (see FIG. 3A), discontinuities (see FIG. 3B), isolations (see FIG. 3C), and the like that are associated with bone loss. [0081] Regions that caused permanent bone loss and/or recoverable bone loss can be identified (116 of FIG. 1). Any suitable technique for identifying permanently lost or recoverable resorption cavities, such as perforations, discontinuities, and isolations, can be utilized. For example, the method described in Saha, P.K. et al., "A new shape preserving parallel thinning algorithm for 3D digital images," Pattern Recognition, Vol. 30, Issue 12, December 1997, ρp.1939-1955, can be utilized.

[0082] The bone formation period and formation deficit can be provided as process parameters (118 of FIG. 1). The bone formation period is a measure of the time for bone formation to occur and the formation deficit is difference between the amount of bone lost from resorption and the amount of bone reformed, expressed as a percentage. These process parameters can be obtained experimentally or from the literature.

[0083] A simulated bone formation process can be carried out in the resorption cavities that did not result in permanent bone loss (120 of FIG. 1), as governed by the parameters specified in 118 of FIG. 1. In carrying out the bone formation process, isolated voxels or voxel-clusters can be removed, as described in Hoshen, J. and Kopelman, R., "Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm," Phys. Rev. B, 1976, vol. 14, pp. 3438-3445. The bone formation process can lead to an updated image of the bone, from which 100 through 120 of FIG. 1 can be repeated, as desired or needed.

[0084] It should be noted that other embodiments are within the scope of the disclosed subject mater. For example, formation may be enhanced near voxels that have the highest stress fields around them, as determined in 106 of FIG. 1. Alternatively, a new FEA can be carried out after 114 of FIG. 1 to determine the location of highest stress fields after the formation of resorption cavities. [0085] In other embodiments, the disclosed subject matter can further be modified to simulate treatment of osteoporosis using clinical drugs, such as bisphosphonate treatment. For example, treatment using bisphosphonate can be simulated by (1) reducing the activation frequency in 108 of FIG. 1, (2) increasing secondary mineralization in the bone trabeculae where the density of individual trabeculae increases, (3) providing a positive bone balance where the thickness of the individual trabeculae can increase over time, sometimes accompanied by decrease in resorption cavity depth.

[0086] Bone morphological analysis (102 of FIG. 1), which will be described in greater detail below, can be carried out at any time after 100 of FIG. 1. Having performed the bone morphological analysis based on the current bone image in 100 of FIG. 1, risk of fracture for the current trabecular bone morphology can be determined (104 of FIG. 1), as will be described in greater detail below. As will be readily apparent to one of ordinary skill in the art, although 102 and 104 of FIG. 1 can be carried out for each iteration of 100 through 120, they need not necessarily be carried out during each iteration. For example, it may be determined that 102 and 104 of FIG. 1 can be carried out after each 2, 5, 10, 15, 20, 50, 100, or any other suitable number of iterations. Moreover, 102 and 104 of FIG. 1 can be carried out at any suitable time during 100 to 120 of FIG. 1 and need not be limited to any specific number of iterations.

[0087] FIG. 4 shows additional details of 102 and 104 of FIG. 1. As shown, morphological analysis of the bone image in 100 can begin with skeletonization of the trabecular bone structure (200 of FIG. 4). Skeletonization can be thought of as an iterative erosion process where bone voxels are peeled off layer-by-layer until no further bone voxels can be removed without altering the shape or topology (e.g., connectivity, tunnels, and cavities) of the trabecular bone architecture. In other words, the trabecular bone image can be reduced to one-dimensional ("ID") or two-dimensional ("2D") structures where the bone may be represented by surfaces and/or curves. FIG. 5 A shows an image of bone obtained using μCT and thresholded as described above in 100 of FIG. 1 and FIG. 5B shows the skeletonized image of the bone in FIG. 5 A. It should be noted that during skeletonization, information regarding resorption cavities, such as perforation, in the trabecular bone can be maintained as shown in FIGS. 6 A and 6B.

[0088] FIG. 7A shows a portion of a bone image obtained in μCT and thresholded as described in 100 of FIG. 1. FIG. 7B shows a skeletonized image of the bone, where the skeletonization is carried out as described in 200 of FIG. 4. As shown in FIG. 7B, individual voxels of the skeletonized image can be classified as curves, surfaces, surface edges, junctions, and the like (202 of FIG. 4). It should be noted that in certain embodiments, the skeletonized image can be further thinned to obtain an arc skeleton where surfaces can be further thinned representatively as arcs. FIG. 7C shows an arc skeleton image where further thinning of surface structures was carried out.

[0089] The classification of voxels into surfaces (or arcs), curves, junctions, edges, and the like can be determined by analyzing the number of components, tunnels, and cavities in a sufficiently small neighborhood of the voxels of interest after each layer removal. To address any potential anomalies that can occur while classification is carried out, the following can be utilized: (1) perform partial classification using the three local topological parameters of components, tunnels, and cavities, (2) classify the voxels into unique classifications by analyzing the results of partial classification over a 3x3x3 neighborhood, and (3) rectify any misclassifications by analyzing the topology of edges. For example, if an arc skeleton image is utilized, and if the length of the arcs and curves in the arc skeleton were measured and found to be shorter than a certain predetermined minimum length criterion, they can be merged with the longest neighboring arc or curve.

[0090] Having classified the voxels in the skeletonized image as curves or surfaces (or arcs), reverse filling of the skeletal structure can be carried out (204 of FIG. 4). For example, just as skeletonization can be thought of as a layer-by-layer erosion of the original structure, reconstruction of the image can be thought of as a layer-by-layer reverse filling of bone regions starting from the set of skeletal voxels until a structure substantially resembling the actual bone structure is obtained (e.g., having a similar bone mass volume, bone density, etc.). During this reconstruction process of reverse filling, the reconstructed voxels can be classified in any suitable manner (206 in FIG. 4). For example, each voxel can be assigned an iteration number where a particular iteration number can denote a particular layer and the iteration number can represent a particular thickness of the layer from the skeletal core. In addition, each voxel can be assigned a topological value using the earlier identified classification scheme. For example, all surface (or arc) related voxels (such as surface edges, surface interiors, surface-surface junctions, etc.) can be assigned a topological value of "1." Similarly, all curve related voxels (such as curve edges, curve interiors, curve-curve junctions, etc.) can be assigned a topological value of "2." For voxels that lie at a junction between a surface and a curve, a topological value of "1.5" can be assigned. As the layers are built up and classified (repeating 204 and 206 of FIG. 4), the topological value of a voxel can be taken as the mean of the topological value of all previously reconstructed voxels in the 3x3x3 neighborhood of the particular voxel. Hence, at the end of the reverse-filling (reconstruction) process (see FIG. 7D), each voxel can have an iteration number based on the number of iterations and a topological value that can range from "1" to "2." As the reconstructed structure is now a 3D structure, voxels having topological value closer to "1" can be classified as more likely to be a plate, and voxels having topological value closer to "2" can be classified as more likely to be a rod. Voxel with a value of 1.5 can be designated to be more like a plate, although they may be designated as rods if more suitable. [0091] It should further be noted that plates and rods can be selectively reconstructed from the trabecular bone skeleton. In plate reconstruction, the voxels labeled as surface can be selectively reconstructed to produce a final structure of plate-type trabeculae that are interconnected by skeletal curves. In rod reconstruction, the voxels labeled as curves can be selectively reconstructed to produce a final structure of rod-type trabeculae that are interconnected by skeletal surfaces.

[0092] FIG. 8A shows an original μCT image (full voxel model) of a trabecular bone sample. FIG. 8B shows a plate-reconstructed image where plates can preserve their thickness or bulk information while rods are represented by skeletal curves. FIG. 8C shows a rod-reconstructed image where the rods can preserve their thickness information while only skeletal structures are preserved for plates.

[0093] A reconstructed image of the bone having voxels classified as plates, rods, and the like has significant advantages in facilitating the prediction of fracture risks based not only on bone mass density, but also on the morphological architecture of the trabecular bone. Hence, the reconstructed image obtained in 206 of FIG. 4 can be utilized to analyze and elucidate important factors that cause fracture in the trabecular bones. For example, morphological and/or finite element analysis of the individual trabeculae and/or bone structure can be performed (208 of FIG. 4), using a variety of reconstructed images to elucidate fracture mechanisms.

[0094] To perform a morphological analysis of individual rod and/or plates of the trabecular bone, the dimensions of the individual rods and plates can be determined as follows. The thickness t of a trabecular plate can be calculated by averaging the local thickness of the trabecula along the surface skeleton. The length t of the trabecular plate can be calculated by £ = ^V_B/t , where VB is the bone volume associated with the trabecula. The length £ of a trabecular rod can be calculated as the length of the curve skeleton. The thickness d of the trabecular rod can be calculated by d = ^4V_B /π£ . For a trabecular bone sample, the average and standard deviation of parameters, such as trabecular plate or rod length, and trabecular plate thickness or rod diameter, can be determined. [0095] A 3D principal component analysis (PCA) technique can be utilized to obtain both geometric and orientation information of trabecular plates and rods (see Sinha, U. and H. Kangarloo, "Principal component analysis for content-based image retrieval," Radiographics, Vol. 22(5), 2002, pp. 1271-1289; Strother, S.C., et al, "Principal component analysis and the scaled subprofile model compared to intersubject averaging and statistical parametric mapping: I. 'Functional connectivity' of the human motor system studied with [15O]water PET," J Cereb Blood Flow Metab, Vol. 15(5), 1995, pp. 738-753; and Tsai, A., et al., "A shape-based approach to the segmentation of medical imagery using level sets," IEEE Trans Med Imaging, Vol. 22(2), 2003, pp. 137-154). The PCA technique is straightforward to apply to bone voxel data.

[0096] To demonstrate, a trabecular plate or rod can be assumed to contain a cluster of voxels with spatial coordinates of Xj, yj, and z; (i=l,M), where M is the number of voxels in this trabecula. First, the center of the trabecula can be calculated as shown in Equations [1] through [3]:

[i]

M

[2]

;=1 M

^Z == L > .z, /M [3]

/=i [0097] Next, a covariance matrix can be constructed as shown in Equation [4]: - Z) z) [4]

[0098] Next, the eigenvalues £_\, I₂, h (^i^>-<?₂ ^>^₃) and eigenvectors H₁, H₂, and H₃ (unit vectors) of this matrix C can be determined. For a trabecular plate, £_\ and H₁ can represent trabecular plate length in the principal direction while ^₃ can represent the plate thickness. For a trabecular rod, £\ can represent the length of the rod while £₂ and ^₃ can represent the rod thickness in the principal directions. [0099] Then, the eigenvalues can be combined with their respective eigenvectors into three vectors, I_x, J₂, and I₃. Therefore, not only can they represent the magnitude of length, but they can also represent the principal direction of individual trabecula. Moreover, these vectors can be averaged to represent not only the average length or thickness of plates or rods, but also to indicate the principal orientation of trabecular plates or rods with the direction of the averaged vector {see FIG. 9). When performing the averaging of these vectors, care may be needed such that the positive direction of the vectors are consistently assigned with respect to one of the sample orientations.

[0100] In certain embodiments, the two length values of trabecular plates as well as the two thickness values of trabecular rods can be further averaged to a single value for a sample. By adopting the American Society of Bone and Mineral Research bone histomorphometric nomenclature {see Parfitt, A.M., et al., "Bone histomorphometry: standardization of nomenclature, symbols, and units. Report of the ASBMR Histomorphometry Nomenclature Committee," J Bone Miner Res, Vol. 2(6), 1987, pp. 595-610), the following morphological parameters can be^" reported: trabecular plate thickness (Tb.pTh), trabecular plate spacing (Tb.pSp), trabecular rod thickness (Tb.rTh), and trabecular rod spacing (Tb.rSp). Trabecular plate or rod numerical density can be calculated by counting the total number of trabecular plates or rods and then dividing by the sample volume. They are named as: trabecular plate number density (Tb.pN/TV) and trabecular rod number density (Tb.rN/TV). The trabecular plate bone volume fraction or trabecular rod bone volume fraction can be calculated by the sum of all the plate or rod volumes divided by the sample bulk volume: trabecular plate bone volume (pBV/TV) and trabecular rod bone volume (rBV/TV). In addition, the rod-rod junction density (Tb.rrN/TV), and rod-plate junction density (Tb.rpN/TV) can be easily calculated by counting the total number of those junctions, respectively, and then dividing by the bulk volume of the sample.

[0101] In addition, anisotropy and anisotropy ratios for trabecular plates and rods can also be calculated: trabecular plate anisotropy l\, i\, l\, and trabecular plate anisotropy ratio DA_p= t\ Ii ₃, trabecular rod anisotropy P₁, P₂, £\, and trabecular rod anisotropy ratio OAf= £\ /P₃.

[0102] Then, multilinear backward regression analyses can be performed to determine the significant morphological determinants among the above-mentioned trabecular bone morphological parameters for the elastic modulus, yield strength, and yield strain of human trabecular bone.

[0103] To perform FEA, in addition to the FEA examples described above, commercially available finite element software, ABAQUS, can be utilized. For each rod trabeculae, two end junctions and a median skeletal point in between the junctions can be utilized to prescribe a 3-node beam element (B32 in ABAQUS), with the capability to also implement curved beam structures {see Guo, X.E. and Kim, C.H., "Mechanical consequence of trabecular bone loss and its treatment: a three-dimensional model simulation," Bone, Vol. 30(2), 2002, pp. 404-411). Such a 3-node beam element can allow for transverse shear deformation, which can be important when a trabecula is short and thick. The beam cross-section of the trabecula can initially be assumed to be circular and uniform along the trabecula. The diameter of the trabecula can be approximated to be equal to the average thickness of the trabecular rod, which can be determined as described above using Equations [1] through [4]. [0104] To perform a triangular meshing for finite element analysis for the plate trabeculae, several points on the plate skeleton can be used. For example, the junctions between surface arcs, the junctions between the curve and the surface, and the surface edge can be used (see FIG. 10). A Delaunay algorithm can be employed for automatic triangular meshing {see Weatherill, N.P., "Delaunay triangulation in computational fluid dynamics," Computers and Mathematics with Applications, Vol. 24(5/6), 2002, p. 129-172). Triangular shell elements can be used to model these resulting triangular plates (S3/3SR in ABAQUS). General purpose shell elements can be used to consider transverse shear deformation, which may be needed for thick trabecular plates. The average thickness of the plate can be used as the thickness of the shell element. In both beam and shell elements, finite strain and large deformation options can be used. The same failure criteria and elastic modulus as in the voxel based models described above can be used. For implementation, a user defined material subroutine (UMAT) can be written for the ABAQUS software. The elastic modulus can be calculated from the initial linear step of the computation and the yield stress and strain can be determined using a 0.2% off-set method. The trabecular bone tissue stress/strain such as the principal stress/strain, strain energy density (SED)₅ and the like can also be extracted for each individual trabecula. The results from the reduced plate/rod models can be compared to those of voxel based models at the level of the overall structure as well as at the individual trabecula level. [0105] Such morphological and/or finite element analyses can elucidate different fracture mechanisms and provide tools for predicting fracture of trabecular bone based on current and simulated/predicted microstructures of trabecular bone.

EXAMPLES

[0106] The following examples provide further detailed description of the disclosed subject matter in accordance with certain embodiments of the disclosed subject matter. The disclosed subject matter can, however, be embodied in many different forms and the following example should not be construed as limiting the scope of the disclosed subject matter in any way.

EXAMPLE 1

[0107] Twenty-nine human trabecular bone samples were obtained from sixteen lumbar vertebrae (75.3 + 13.8 year-old, 8 males and 8 females), four proximal femurs (one pair: 60 year-old female, two singles: 64 year-old and 44 year-old males) and one pair of proximal tibiae (69 year-old male). The subjects were screened to exclude metabolic bone diseases or bone cancer, and X-ray radiographs were taken to ensure that there was no evidence of damage or other bone pathologies. Twenty 8 mm diameter on-axis (along the principal trabecular orientation) vertebral specimens along the superior-inferior direction were obtained following previously published protocols (see Kopperdahl, D.L., and Keaveny, T.M. "Yield strain behavior of trabecular bone," J Biomech Vol. 31, 1998, pp.601-608). Four on- axis specimens from human proximal tibiae and five on-axis specimens from proximal femurs were obtained using similar protocols (see Chang, W.C., Christensen, T.M., Pinilla, T.P., Keaveny, T.M. "Uniaxial yield strains for bovine trabecular bone are isotropic and asymmetric," J Orthop Res., Vol. 17, 1999, pp. 582-585; Morgan, E.F., Keaveny, T.M. "Dependence of yield strain of human trabecular bone on anatomic site," J Biomech Vol. 34, 2001, pp. 569-577; and Keaveny, T.M., Guo, X.E., Wachtel, E.F., McMahon, T.A., Hayes, W.C. "Trabecular bone exhibits fully linear elastic behavior and yields at low strains," J Biomech, Vol. 27, 1994, pp. 1127-1136). After thawing at room temperature, the cylindrical specimens were aligned and stabilized with wet gauzes in a 15 ml centrifuge tube along their longitudinal axis. The specimen tube was fastened in the specimen holder of a μCT system, vivaCT 40 (SCANCO Medical AG, Bassersdorf, Switzerland). The central gauge length of 15 mm was scanned at 21μm nominal isotropic resolution and the central ~4x4x4 mm cubical sub-volume equivalent to about 191x191x191 voxels was extracted from each reconstructed image. A global thresholding technique was applied to binarize gray-scale μCT images where the minimum between the bone and bone marrow peaks in the voxel gray value histogram was chosen as the threshold value. Isolated voxels or disconnected voxel-clusters were removed from binarized μCT images by computing the largest bone component under 26-connectivity (see Saha, P.K., Chaudhuri, B.B. "Detection of 3-D simple points for topology preserving," IEEE Trans Pattern Anal Machine Intell, Vol. 16, 1994, pp. 1028-1032; and Hoshen, J., Kopelman, R. "Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm," Phys Rev B, Vol. 14, 1976, pp. 3438-3445). The resulting μCT images (see, e.g., FIG. HA) were used for digital topological analysis (DTA) and FEA.

[0108] Thereafter, skeletonization (200 of FIG. 4) of the μCT images was carried out as shown in FIG. HB. Digitial topological classification was carried out (202 of FIG. 4) where each skeletal voxel was uniquely classified as a surface, curve type, edge, junction, or the like. As shown in FIG. HC, the surface voxels are shown as dark grey whereas the curved voxels are shown in lighter shading. Next, reconstruction of the skeletonized image (204 and 206 of FIG. 4) was carried out. FIGS. HD through HF show various types of reconstructed models. For example, FIG. HD shows a fully reconstructed structure with the plate voxels shown as dark gray and rod voxels shown in lighter shading. FIG. HE shows a rod- reconstructed model where the surface voxels are shown in dark gray and the rod voxels are shown in lighter shading. FIG. HF shows a plate-reconstructed model where the plate voxels are shown in dark gray and the curve voxels are shown in lighter shading. Skeletonization, DTA, and reconstruction were performed using a software written in MICROSOFT VISUAL C⁺⁺ (MICROSOFT FNC, Redmond, WA, USA) running on a DELL XPS PC workstation (DELL INC., Round Rock, TX, USA).

[0109] Four different μCT images of each specimen: original, skeletonized, rod- reconstructed, and plate-reconstructed images were converted to a set of voxel-based FE models by converting each voxel to an 8-node brick element. The original, full voxel model had sufficient spatial resolution to appropriately characterize elastic modulus of trabecular bone. The other three models were synthetic models with specific purpose of examining contributions of micro-architecture or trabecular type (rod vs. plate). For each model, a linear FE analysis was applied to determine the apparent Young's modulus E^* using an element-by- element pre-conditioned conjugate gradient solver (see Hollister, SJ., Brennan, J.M., Kikuchi, N. "A homogenization sampling procedure for calculating trabecular bone effective stiffness and tissue level stress," J Biomech, Vol. 27, 1994, pp. 433-444).

[0110] The trabecular bone tissue was modeled as an isotropic, linear elastic material with a Young's modulus E_s of 15 GPa and a Poisson's ratio of 0.3 {see Guo, X.E., Goldstein, S. A. "Is trabecular bone tissue different from cortical bone tissue?" Forma, Vol. 12, 1997, pp. 185-196). Apparent modulus values were determined along three orthogonal directions: the principal trabecular orientation (z) and other transverse directions (x and y). In each direction, a fixed

mm (i=x, y, and z) was applied perpendicularly to one face of the model while the opposite face was imposed zero displacement along the same direction. The total reaction force (RFi) was calculated from the FE analysis and the apparent Young's

Ti Tj n modulus E^* by: E- = *-, where t_λ and A; were the length and cross-sectional area of the

specimen in the i^th direction (i=x, y, and z).

[0111] To further describe the microstructure of trabecular bone, the following morphological parameters were also calculated from the labeled bone images: bone volume fraction (BV/TV) (the total volume of bone voxels divided by the bulk volume); plate bone volume fraction (pBV/TV) (the total volume of plate bone voxels divided by the bull- volume); rod bone volume fraction (rBV/TV) (the total volume of rod bone voxels divided by the bulk volume); the plate-to-rod ratio (pBV/rBV) (the total volume of plate bone voxels divided by the total volume of rod voxels); plate tissue fraction (pBV/BV) (the total volume of plate bone voxels divided by the total volume of bone voxels); rod tissue fraction (rBV/BV) (the total volume of rod bone voxels divided by the total volume of bone voxels); and the thinning ratio (the total volume of bone voxels divided by the total volume of skeleton voxels).

[0112] All the statistical analyses were performed using the STATISTIC TOOLBOX in MATLAB 6.0 (THE MATFfWORKS, Inc., Natick, MA, USA) on a PC Workstation. To evaluate the contribution of the micro-architecture to the mechanical properties of trabecular bone, linear correlations between the Young's modulus (E*) under compression of full voxel model and that of the skeleton model were separately studied along all three image coordinate directions. Second, to evaluate the effect of trabecular type on the elastic modulus of entire trabecular bone network, linear correlations between Young's modulus of the full voxel model as well as those of each of the rod- and plate-reconstructed models were separately examined along the three image coordinate axes. The correlation coefficients between BV /TV of foil voxel model and those of each of the skeleton, rod-reconstructed, and plate- reconstructed models were also determined. For the foil voxel model, correlations between Young's modulus along principal orientation of trabecular bone specimen (E^*) and BV/TV as well as pBV/TV and rBV/TV, which represent the trabecular plate and rod density, were computed. Associations between E^* of the foil voxel model and trabecular plate tissue fraction (pBV/BV), rod tissue fraction (rBV/BV), and plate-to-rod ratio (pBV/rBV) were evaluated as well. Finally, correlations between thinning ratio and BV/TV as well as E^* of the foil voxel model were studied.

[0113] To differentiate the role of trabecular architecture from bone volume fraction, partial correlation analyses between elastic modulus of foil voxel model and those of each of the skeleton, rod-reconstructed, and plate-reconstructed models along three coordinate axes with the effect of BV/TV removed were performed. Similarly, partial correlation between plate (rod) tissue fraction and E^* of the foil voxel model with effect of BV/TV removed was studied.

[0114] As shown in FIG. 12, the plots for the skeleton model and rod-reconstructed model were almost superimposed on each other. Moreover, although the minimal bone volume was maintained in the skeletonized model, surprisingly there was a significant correlation between E^* of the skeleton model and that of the foil voxel model along the principal orientation of trabecular bone specimens. Significant correlations were also found in the other two transverse directions with slightly different correlation coefficients as shown Table 1. Although, the elastic moduli of the skeleton models correlated significantly with those of the original, foil voxel models, the values of the skeleton moduli were more than one order of magnitude less than the foil moduli values. This, of course, was because the skeletonization led to, on average, a factor of 8.91 reduction in bone volume fraction. Moreover, partial correlation analysis revealed that the correlations between elastic modulus of foil voxel model and those of each of the skeleton, rod-reconstructed and plate- reconstructed models were highly significant even with the effects of BV/TV removed (Table 1). Table 1: Correlation coefficients between Young's modulus of trabecular bone

: p<0.005 for all correlations in Table 1. ⁺: partial correlation with the effect of BV/TV removed.

[0115] As described above and shown in FIG. 12, when trabecular rods were fully reconstructed to their original bone volume in the rod-reconstructed models, the changes in E_z ^* from the skeleton model were negligible in terms of the correlation between E^* of the rod-reconstructed model and that of full voxel model. On the average, the rod-reconstructed model had a bone tissue volume, which corresponded to 21.3% of the original bone tissue volume. But the rod-reconstructed model only represented 1.5% of the original elastic modulus. The results in the other two transverse directions were similar (see Table 1) [0116] On the other hand, FIG. 12 shows that when trabecular plates were fully reconstructed to their original bone volume in the plate-reconstructed models, the changes in elastic modulus from the skeleton model were quite dramatic. The plate-reconstructed models had 90.3% of the original bone tissue volume and 53.2% of the original elastic modulus on average. The results in the other two transverse directions were similar (see Table 1).

[0117] These observations in elastic modulus are consistent with the results in bone volume fraction, shown in FIG. 13, which suggest that the plate-reconstructed model maintained most of the bone volume fraction compared with the rod-reconstructed and skeletonized model.

[0118] The trabecular bone microstructure in the specimens studied ranges from rod- dominant to plate-dominant with the plate-to-rod ratios ranging from 4.6 to 29.7. FIGS. 14A through 14C shows the full voxel models and skeleton models of three trabecular bone samples, illustrating the variations in structure type from plate-dominant to rod-dominant. The left of each figure shows original μCT image (full voxel model) of the trabecular bone sample while the right shows skeletonized image (skeleton model). FIG. 14A has a plate-to- rod ratio of 20.10, plate tissue fraction of 0.9526, and a Young's modulus of 974.8 MPa. FIG. 14B has a plate-to-rod ratio of 10.22, plate tissue fraction of 0.9109, and a Young's modulus of 674.8 MPa. FIG. 14C has a plate-to-rod ratio of 5.94, plate tissue fraction of

0.8559, and a Young's modulus of 216.8 MPa. [0119] The elastic moduli E^* of these samples were also determined and indicated a direct relationship with the plate-to-rod ratio. In FIG. 15A, E^* of the full voxel model correlated significantly with the plate bone volume fraction (pBV/TV) by a nonlinear regression of power law while no significant correlation was found between E^* and rod bone volume fraction (rBV/TV) (FIG. 15B). It should be noted that the pBV/TV combines trabecular type and bone volume fraction. To isolate the effect of the bone volume fraction, pBV/TV and rBV/TV were normalized by BV/TV to obtain plate tissue fraction (pBV/BV) and rod tissue fraction (rBV/BV). Plate or rod tissue fraction represents the percentage of bone tissue belonging to the trabecular plate or rod and the sum of them always equals 1. The correlation between pBV/B V and E^* , and between rBV/BV and E_z ^* are shown in FIG. 15C and FIG. 15D to follow nonlinear regression of power laws. Both correlations were significant (r²=0.49 & 0.42, p<0.001) but have opposite sign. Furthermore, partial correlation analysis demonstrated plate fraction significantly correlated to Young's modulus of full model (r²=0.22, p=0.01) with the effect of BV/TV removed. Finally, correlations between E] and BV/TV, and between E] and plate-to-rod ratio are shown in FIG. 15Ε and FIG. 15F. [0120] The partial correlation between the elastic moduli of skeleton models and those of full models suggests that the skeleton of the trabecular bone, which is the representative of the network's topology, affects its mechanical properties in a manner independent of bone volume fraction. Certainly, bone volume fraction is also a significant factor, which determines elastic modulus of trabecular bone (FIG. 15E) as bone volume fraction is also shown to be highly correlated with the elastic modulus. However, FIG. 13 shows that although the skeleton has lost most of the bone volume, BV/TV of skeleton model is still significantly correlated with full BV/TV. Thus, the data suggests that the micro-architecture of trabecular bone may have a direct relationship to BV/TV, which is consistent with previous 2D analyses {see Chappard, D., Legrand, E., Haettich, B., Chales, G., Auvinet, B., Eschard, J.P., Hamelin, J.P., Basle, M.F., Audran, M., "Fractal dimension of trabecular bone: comparison of three histomorphometric computed techniques for measuring the architectural two-dimensional complexity," J Pathol., Vol. 195, 2001, pp. 515-521; and Thomsen, J.S., Ebbesen, E.N., Mosekilde, L., "Relationships between static histomorphometry and bone strength measurements in human iliac crest bone biopsies," Bone, Vol. 22, 1998, pp. 153- 163). It should be noted that the strong correlation of the elastic moduli with both the bone volume fraction and the microarchitecture suggests that the bone volume fraction may depend on the microarchitecture of the bone.

[0121] Moreover, the plate-reconstructed model represents a hypothetical trabecular bone structure, which has lost most of the rod-like trabecular volume information while maintaining all plate-like trabeculae and the skeleton of rod-like trabeculae. Contrarily, the rod-reconstructed model keeps all rod-like trabeculae while preserving the skeleton of plate- like trabeculae. These results further show that the loss of trabecular rod volume does not significantly affect either elastic modulus or overall bone volume fraction compared with the loss of plate volume (see FIGS. 12 and 13). Moreover, the elastic modulus is significantly correlated with pBV/TV while no correlation exists between the elastic modulus and rBV/TV (FIGS. 15A and 15B). Hence, the trabecular plates may play a much more important role in determining the elastic modulus and bone volume fraction than do rods. [0122] Furthermore, plate and rod tissue fractions quantitatively indicate the percentage of the bone tissue that belongs to either trabecular type. In this manner, the quantitative influence of trabecular types on elastic modulus of trabecular bone can be demonstrated (FIGS. 15C and 15D). An alternative approach toward characterizing trabecular bone type is the structural model index (SMI), which is a measure of the bone's degree of plate-likeness or rod-likeness in the form of an index ranging from 0 to 3 {see Hildebrand, T., Ruegsegger, P., "A new method for the model independent assessment of thickness in three-dimensional images," J Microscopy, Vol. 185, 1997, pp. 67-75). The SMI for all 29 samples were evaluated using the software on the vivaCT 40 (SCANCO Medical AG, Bassersdorf, Switzerland) and they were compared with the plate-to-rod ratio. A significant linear correlation existed between the SMI and the plate-to-rod ratio (r² = 0.85, p<0.001), further corroborating the importance of plate- type trabeculae. It should be emphasized that the plate- to-rod ratio has an advantage of identifying structural elements as belonging to a particular trabecular type rather than only a statistical calculation of SMI.

EXAMPLE 2

[0123] A set of fifteen human trabecular bone samples consisting of six vertebral, five femoral, and four tibial cores (~8 mm in diameter and 30 mm in length) were scanned at 21 μm resolution using a vivaCT 40 μCT system. The central ~4x4x4 mm cubical sub-volume was selected from each reconstructed μCT image and a global thresholding followed by a clustering analysis was applied to obtain the trabecular bone region shown in FIGS. 16A and 16F within each sub-volume.

[0124] A morphological analysis was carried out as follows. First, a skeletonization and topological classification algorithm based on digital topological analysis (DTA) was applied to first iteratively peel bone surface layers until only the medial surface skeletons, or curve skeletons, are left (200 of FIG. 4). Then, a local topological analysis was performed to determine the class of each skeletal voxel (202 of FIG. 4). FIGS. 16B and 16G show the skeletal images of the trabecular bone. Each surface skeleton was further thinned to an arc skeleton using DTA-based arc-thinning algorithm (FIGS. 16C and 16H). The branching points were then identified in the arc-skeleton and were used to split the skeleton into individual arcs. Spurious branches determined in the initial arc-skeleton based on their lengths were merged with the longest neighboring arc or curve. Final partitions of individual trabeculae were computed by iteratively growing these individual arcs within the trabecular region while retaining the topological class of each voxel (204 and 206 of FIG. 4). Each such partition was treated as a trabecula (FIGS. 16E and 16J), and the normal of each plate- type trabecula was determined by fitting a plane to its surface skeleton and the thickness was measured along the normal. The orientation of each rod-type trabeculae was obtained by fitting a line to its curve skeleton and the thickenss was calculated from its volume and length. [0125] pTb.Th (mm) and rTb.Th (mm) represented average plate and rod thickness. Total plate trabeculae number and rod trabeculae number was normalized by the sample volume to get the pTb.N (I/mm) and rTb.N (I/mm). Plate bone volume (pBV) and rod bone volume (rBV) was calculated and normalized by the total volume (TV) to obtain the relative plate fraction (pBV/TV) and rod fraction (rBV/TV). Junction number of skeletonized image was normalized by sample volume to obtain the junction density (I/mm). [0126] μCT images were also converted to voxel-based finite element models and subjected to linear analyses to determine the Young's modulus (E). The trabecular bone tissue was modeled as an isotropic, linear elastic material with a Young's modulus E_s of 18 GPa and a Poisson's ratio of 0.3. All the new structural indices were compared with standard 3D structural indices by linear regression and all the new and standard indices were correlated to the elastic moduli and BVF using a linear regression. [0127] The BVF of 15 trabecular bone samples varies from 8.6% to 18.3% with the Young's modulus ranging from 260 to 1170 MPa. Table 2 shows the basic statistics of the indices pTb.Th, rTb.Th, pTb.N, rTb.N, pBV/TV, and junction density. Table 3 shows the basic statistics of standard 3D structural indices trabecular number (Tb.Th), trabecular thickness (Tb .N), structure model index (SMI) and connectivity {see Kabel, J. et al,

"Connectivity and the elastic properties of cancellous bone" Bone, 1999, Vol. 24(2), pp. 115-

20).

Table 2: Structural indices based on individual trabeculae segmentation.

Table 3: Standard structural indices based on individual trabeculae segmentation.

[0128] Significant correlation by linear regression was found within four groups: (1) pTb.Th and Tb.Th, (2) pTb.N and Tb.N, (3) pBV/TV and SMI, (4) junction density and connectivity (Table 4). No correlation is found between rTb.Th and Tb.Th. as well as between Tb.N and rTb.N.

[0129] pTb.Th, pTb.N and pBV/TV have relatively weak or no correlation with BVF compared with Tb.Th, Tb.N and SMI.

[0130] Two groups: (1) Tb.N and pTb.N and (2) SMI and pBV/TV correlated with elastic moduli.

Table 4: Correlation coefficient (R²) between new and standard structural indices and between bone volume fraction, elastic modulus and all the structural indices. *p<0.001

[0131] The morphological analysis technique described above allows segmentation of trabecular bone μCT images into the individual trabeculae pieces as trabecular plates or rods. Therefore the standard morphological features such as Tb.Th and Tb.N can be further described as pTb.Th, rTb.Th, pTb.N and rTb.N.

[0132] As shown above, pTb.Th and pTb.N highly correlates with Tb.Th and Tb.N while rTb.Th and rTb.N shows no correlation. This implies the plate-type trabeculae dominate in the architecture of 15 samples, which is further confirmed by an averaged 0.9 plate fraction (pBV/TV).

[0133] The index pBV/TV, which can quantitatively evaluate the plate volume ratio, shows a weak correlation with BVF but a significant correlation with elastic modulus. This result not only confirms that a lower bone mass is more likely characterized by more rod-like structure and also suggests that other than bone volume fraction the volume distribution of plate-like and rod-like structure also plays important role in determining mechanical properties of trabecular bone.

EXAMPLE 3

[0134] Fifteen cylindrical (~8 mm in diameter and 20 mm in length) human trabecular bone cores were cored from L4 (n=14) and L5 (n=l) vertebral bodies of fifteen donors (6 males/9 females; age 74±11). 3D high resolution images were obtained for each specimen using μCT (Scanco μCT 20; Scanco Medical AG, Switzerland). The central ~5x5x5 mm cubical sub-volume was extracted from each reconstructed image and converted into a voxel- based finite element (FE) model of 44 μm resolution. Geometrically non-linear FE analysis using the nonlinear material properties having tension-compression strength asymmetry of the hard tissue, were conducted on an IBM SP4 parallel supercomputer (IBM Corporation, Armonk, NY) to perform ten steps evenly increasing compression from 0 up to 1% apparent strain. A finite deformation based plasticity model was used to model trabecular tissue with previously calibrated yield properties and a post-yield modulus equal to 5% of the initial modulus {see Bevill, G. et al, "Large deformation effects in failure behavior of trabecular bone" ORS, 2005; and Bayraktar et al, "Mechanisms of uniformity of yield strains for trabecular bone", J. Biomech., 2004, Vol. 37, pp. 1671-78). A custom parallel FE solver was used to solve these problems in parallel {see Admas et al, "Ultrascalable implicit finite element analyses in solid mechanics with over a half a billion degrees of freedom" ACM/IEEE Procs. SC2004, 2004). The voxel element, which reached yielding, was marked at each loading step. [0135] A digital topological analysis (DTA) {see Saha, P. K. et al, "A new shape preserving parallel thinning algorithm for 3D digital images" Pattern Recogn., 1997, VoI 30(12), pp. 1939-55) based segmentation technique was applied to segment the trabecular bone structure into individual trabeculae following their topology. The trabecular type, volume and orientation (horizontal (0..θ ≤$0°), oblique (30°<θ ^>0°) or longitudinal (60°<θ -_90°), where θ is the angle with respect to horizontal direction) were identified for each segmented plate trabecula (FIG. 17A) and rod trabecula (FIG. 17B). Yielding of any voxel within a trabecula was defined as the point of failure. To get the percentage of yielded trabecular plate and rod, at each loading step, the number of yielded trabecular plates and rods was normalized by the total number of each type trabeculae in the sample (FIGS. 18A and 18B), and these data were stratified by orientation (FIGS. 19A and 19B). The percentage of yield volume in each type of trabecula was also quantified at each step for both trabecular plate and rod (FIG. 18B). Paired t-tests were performed with p ^).O5 considered significant.

[0136] The average bone volume fraction of fifteen samples was 0.105±0.036. Despite the low average bone volume fraction, there were significantly more plates than rods while the majority of plates was in the longitudinal direction and most rods were in the horizontal direction (Table 5). Table 5: Summary of trabecular plate/rod classification

[0137] As shown in FIG. 18, trabecular bone yielding initiated in trabecular rod and then more plates yielded at the later stage (>0.6% global strain).

[0138] As shown in FIG. 19A, among the yielded plates, significantly more longitudinal plates yielded as compared to horizontal and oblique plates through the entire failure process.

[0139] As shown in FIG. 19B, among the yielded rods, significantly more horizontal rods than longitudinal and oblique rods yielded above about 0.3% global strain.

[0140] These results suggest that trabecular bone yield may be initiated in trabecular rods while more trabecular plates yield in the late stage of trabecular bone failure process. This would suggest that the weakest link in trabecular bone is trabecular rods while trabecular plates dominate the overall strength. The results also show that more trabecular rods yield in the horizontal direction than in the other two directions, demonstrating the importance of horizontal rods in the strength of trabecular bone. On the contrary, more trabecular plates yield in longitudinal direction than in the other two directions, demonstrating important roles of longitudinal plate.

EXAMPLE 4

[0141] This example quantitatively studied the relative contribution of trabecular rods in different orientations to elastic moduli of vertebral trabecular bone.

[0142] Fourteen cylindrical (~8 mm in diameter and 20 mm in length) human trabecular bone cores in the axial direction were harvested from L4 (n=13) and L5 (n=l) vertebral bodies of fifteen donors (6 males/9 females; age 74±11). Samples were scanned at 21μm resolution using a micro-computed tomography (μCT) system (μCT 40, Scanco Medical Inc.). The scanned image was thresholded, skeletonized, classified, and reconstructed with both plates and rods as described above. The orientation of the rod was defined to be vertical (0 ≤S <30°; FIG. 20A), oblique (30°<θ ≤S0°; FIG. 20B), or horizontal (0 ≤θ ≤30°; FIG. 20C) according to the angle θ between the macroscopic axial direction and the rod. [0143] Next, vertical rod trabeculae, oblique rod trabeculae, or horizontal rod trabeculae were artificially and selectively removed from the original 3D μCT images of trabecular bone. Four sets of specimen-specific μCT based micro finite element (μFE) models were constructed for all specimens: (1) original; (2) vertical rods removed; (3) oblique rods removed; and (4) horizontal rods removed. In all images and models, trabecular plates were untouched. The bone tissue properties were chosen as isotropic, linear elastic with a Young's modulus of 15 GPa and a Poisson's ratio of 0.3 for all models.

[0144] Using an element-by element pre-condition conjugate gradient solver, six μFE analyses were performed for each model, representing three compression tests and three shear tests, resulting in a total of 336 FE analyses. The full anisotropic stiffness tensor was calculated first and then elastic material constants (three Young's moduli, Exx, Eyy, Ezz and three shear moduli, Gyz, Gxz, Gxy,) were derived (B. Van Rietbergen, et al., J. Biomech. 1996, VoI 29(12), pp. 1653-57). BV/TV, plate volume fraction (pBV/TV), rod bone volume fraction (rBV/TV), and rod number (rod #), were also calculated for all models. Student paired t-tests were performed between four types of models with a significant level of p < 0.05.

[0145] Average BVF, rBV/TV, and rod # of fourteen specimens were calculated for each category of models based on a morphological analysis (Table 6). Trabecular rods accounted for 15.1% of the total BV/TV while the pBV/TV remained constant of 0.09 (84.9%) for all four models, as trabecular plates were unchanged. Horizontal rods represent 62.5% of all rods in terms of bone volume, and 58.2% of all rods in terms of rod number. Vertical rods represent 12.5% of all rods in terms of bone volume, and 12.2% of all rods in terms of rod number. Oblique rods represent 31.3% of all rods in terms of bone volume, and 29.7% of all rods in terms of rod number.

Table 6: Average BVF, rBV/TV, and rod # for original model and 3 reduced models.

[0146] As shown in FIG. 21, after removing the vertical rods, oblique rods or horizontal rods, all of the material constants of reduced models decreased significantly compared to the full voxel model. Removing horizontal rods resulted in the most reductions in elastic moduli especially in Exx, Eyy, and Gxy. Removing vertical rods results in the least reductions in elastic properties with the largest reduction in Ezz. Oblique rods seem to contribute equally to all elastic moduli.

EXAMPLE 5

[0147] In this example, a rigorous 3D digital topological analysis (DTA) technique was incorporated in the 3D simulation of trabecular bone remodeling with accurate considerations of three types of microscopic bone loss in the pathophysiology of osteoporosis: perforation

(FIG. 22A), breakage (FIG. 22B), and disconnection (FIG. 22C) of trabeculae.

[0148] Furthermore, a dynamic bone remodeling physiological process simulation, with all parameters derived from recent published clinical data (see Recker, R., et al, "Bone remodeling increases substantially in the years after menopause and remains increased in older osteoporosis patients," JBMR, Vol. 19(10), 2004, p. 1628), was carried out in 3D human trabecular bone models. To explore the contributions of different bone loss mechanisms, the amount of the bone loss due to each mechanism was quantified in the case of menopausal osteoporosis. In addition, the changes in architecture and elastic modulus of trabecular bone during osteoporosis were quantified.

[0149] Twelve human trabecular specimens were obtained from a third lumbar vertebra (52 y.o. male), three proximal femurs (three singles: 64 y.o. male, 44 y.o. male, and 60 y.o. female) and one proximal tibia (69 y.o. male) and scanned at 21μm resolution using a vivaCT 40 μCT system. The central ~4x4x4 mm³ cubical sub-volume was selected from each reconstructed μCT image and globally thresholded.

[0150] A bone remodeling cycle corresponded to a 200-day period in real life (40 days resorption/160 days formation). Resorption cavities (42 μm deep and 126 μm in diameter) were created in 40 days according to the current activation frequency and distributed randomly over the bone surface. Each cavity went through a digital topological detection and the cavity that caused a perforation, breakage, or disconnection of a trabecula was identified. Every resorption cavity was refilled in 160 days unless it caused a perforation or a breakage of a trabecula. Disconnected trabeculae were removed from the simulation. In the meantime, new resorption cavities were continuously created in ,40 days. The simulation was started 5 years before and ended 15 years after menopause. As shown in FIG. 23, the bone remodeling activation frequency was taken to linearly rise from 0.13/year at 5 year before menopause to 0.24/year at 2 years after menopause, to increase to 0.37/year at 7 years after menopause, and then to remain elevated at the same level.

[0151] For each sample, 3D images of -5 (FIG. 24A), 0 (FIG. 24B), 5 (FIG. 24C), 10 (FIG. 24D), and 15 years of remodeling simulation were converted to voxel-based finite element models and subjected to linear analysis to determine the Young's modulus (E) along the principal trabecular orientation. Based on a digital topological analysis, the plate fraction (pBV/TV) was also calculated as described above for those images. [0152] The bone volume fraction (BVF) of eight trabecular bone samples varied from 10.5% to 34.9% with the Young's modulus ranging from 0.337 to 3.31 GPa. FIGS. 24A through 24D show the typical simulation of bone loss during menopause. [0153] The time courses of averaged, relative BVF due to various types of bone loss and the averaged, relative elastic moduli are shown in FIG. 25. The results suggest that the primary bone loss was due to the perforation of trabecular plate and secondarily due to broken trabeculae. [0154] As the activation frequency was increased for each iterative cycle (see FIG. 16), the trabecular bone remodeling simulation predicted decreased BVF, E, and pBV/TV as shown in Table 7 below, which is consistent with previous clinical study, such as that described by Recker, R. et al. in the Journal of Bone and Mineral Research (2004) cited above.

Table 7: The average bone volume fraction (BVF), Young's modulus (E), and plate fraction (pf) during 20 years of perimenopausal remodeling simulation.

[0155] The averaged BVF of samples from two anatomic sites, femoral neck (FIG. 26A) and spine (FIG. 26B), were evaluated separately and compared with previous clinical data. For each group, the pattern of predicted bone loss was consistent with its corresponding clinical data.

[0156] Table 8 shows three representative samples of different initial BVF. Clearly, trabecular bone with higher BVF and pBV/TV underwent a slower bone loss, given same activation rate and menopausal duration.

Table 8: Baseline (5 years before menopause) and predicted relative reduction of BVF, Pf and elastic modulus at 15 years after menopause for three bone samples from femoral neck, tibia and spine.

[0157] As shown above, with all simulation parameters based on physiological measurements, the simulation predicted a bone loss time-course which is consistent with clinical observations. The results of this study suggest that the trabecular plate perforation accounts for more than 70% of bone loss during menopause. Furthermore, the decrease of plate fraction during simulation quantitatively confirmed that transition of trabecular plate into rod is involved with post-menopausal osteoporosis. Bones with higher initial BVF show a slower bone loss time course which is consistent with the fact that peak bone mass is one of the important risk factors in osteoporosis.

[0158] In addition, different bone samples experienced different BVF, pBV/TV, and E reduction rate, which implies that the trabecular bone deterioration due to remodeling is dependent on the initial stage of the bone mass and architecture. [0159] Activation frequency also appears to be a main factor dominating bone remodeling rate. Following the increase of the activation frequency during and after menopause, FIG. 25 suggests that more severe bone loss starts about 5 years after menopause. [0160] Hence, these results provide additional insights into the mechanism of osteoporosis and provide a 3D prediction tool for individuals associated with risk of osteoporosis before they occur.

EXAMPLE 6

[0161] In this example, the contributions of bisphosphonate treatment was evaluated using a three dimensional (3D) simulation of bone remodeling. Bisphosphonate is an antiresorptive agent and has been widely used to treat osteoporosis. Several clinical trials showed that the daily oral bisphosphonate induced a sustained decrease in biochemical marker of bone turnover and increased bone mineral density (BMD) at most skeletal sites in postmenopausal women {see Devogelaer, et al. Oral alendronate induces progressive increases in bone mass of the spine, hip, and total body over 3 years in postmenopausal women with osteoporosis. Bone, 1996, Vol. 18, pp.141-50; and Bone, et al, "Ten years' experience with alendronate for osteoporosis in postmenopausal women", N EnglJ Med, 2004, Vol. 350, pp. 1189-99), which is significantly associated with reduction in vertebral and nonvertebral fractures {see Black, et al, "Fracture risk reduction with alendronate in women with osteoporosis: the fracture intervention trial", J Clin Endocrinol Metab 2000, Vol. 85, pp. 4118-24).

[0162] At tissue level, the action of bisphosphonate appears to be through three different mechanisms: (1) reduced remodeling activation frequency (Ac.f); (2) increased secondary mineralization (MIN); and (3) positive bone balance (PBB) by reduced depth of resorption cavity and increased thickness of newly formed bone. However, the relative contributions of each mechanism have not been quantified at the microstructural level of trabecular bone. [0163] Eight human vertebral trabecular samples were obtained from lumbar vertebrae of seven donors (1 female/7 males, age 68±16) and scanned at 21 μm resolution using a micro computed tomography (μCT) system (vivaCT 40, Scanco Medical AG, Switzerland). The center ~4x4x4 mm cubical sub-volume from each specimen was extracted and thresholded. [0164] The simulation included a menopausal phase followed by a treatment phase. A bone remodeling cycle included a 40-day bone resorption period followed by a 160-day bone formation period. Resorption cavity was 126 μm in diameter and 42 μm in depth and randomly located over the bone surface according to the current activation frequency. When a resorption cavity caused any perforation, breakage, or disconnection of trabeculae, it was assumed new bone does not form to fill the cavity (see Liu, et al, "Simulating 3D architectural and mechanical changes in human trabecular bone during menopause", Trans. 52^nd ORS, 2006). Before simulation, each bone voxel was initialized to 100% MIN. While menopausal phase started, each newly formed bone voxel was reassigned a 70% MIN and the t secondary MIN increased by ( 1- e ^λ ), λ=2.2 year. During the menopausal phase, each bone sample underwent 6 years of bone remodeling (Ac.f = 0.4/year (see Recker, et al, "Bone remodeling increases substantially in the years after menopause and remains increased in older osteoporosis patients", JBMR, 2004, VoI 19(10), pp. 1628-33)) where a realistic MIN distribution in trabecular tissue was achieved. The treatment phase started immediately after the menopausal phase and continued for 10 years. Three mechanisms of bisphosphonate treatment were considered. First, Ac.f was decreased by 90% to 0.04/year (reduction of remodeling space, RS). Second, the increased MIN with time was considered. Finally, a PBB was considered by decreasing the depth of resorption cavity to 21 μm. For each sample, the simulated images of normal (beginning of menopausal phase), menopausal (end of menopausal phase) and treated (end of treatment phase) bones went through a morphological analyses based on digital topological analysis (DTA) (see Saha, et al, "A new shape preserving parallel thinning algorithm for 3D digital images", Pattern Recogn., 1997, Vol. 30(12), pp. 1939-55). Averaged MIN of untreated and treated bones were also calculated and compared.

[0165] The time course of BMD change attributed to different mechanisms is shown in FIG. 27A. As shown, simulations predicted corresponding clinical data when all three mechanisms of bisphosphonate: reduction in remodeling space (RS), increased mineralization (MIN), and positive bone balance (PBB) were considered. Comparing with baseline, untreated menopausal samples from the simulation had a significant decrease in bone volume fraction (BV/TV), BMD, plate fraction, trabecular plate and rod thickness (pTb.Th and rTb.Th) and significant increase in trabecular plate and rod number density (pTb.N and rTb.N) (Table 9).

Table 9: The averaged BV/TV, BMD change, MIN and morphological parameters of normal, menopausal and treated bone specimens.

[0166] During ten years of simulation of the bisphosphonate treatment, BV/TV of trabecular bone was restored to a similar level as before the menopause. As shown in Table 9, this increase in BV/TV was achieved mostly via significant increases in pTb.Th, and rTb.Th. The increase in pTb.N during menopause was completely restored by treatment while rTb.N continued to increase. FIGS. 27B and 27C show a 2D distribution of secondary mineralization in trabecular bone before and after treatment, respectively. As shown, bisphosphonate treated bone showed a clear shift toward a higher mineralization level (see also Table 9). Changes in both microstructure and mineralization resulted in the increase in BMD as shown in FIG. 27A.

[0167] As shown, the simulation predicts fairly well the increased BMD following bisphosphonate treatment. Therefore, the results from this study suggest that all three mechanisms (reduced remodeling activation, increased mineralization, and positive bone balance) are important. Notably, this 3D simulation suggests that a positive bone balance is necessary for a long term gain in BMD (FIG. 27A, solid curve). This simulation can also model the mineralization process in bone remodeling and quantitatively predict the increased mineralization during bisphosphonate treatment, which contributes significantly to the increased BMD. [0168] As described above, the disclosed subject matter permits computationally efficient methods for simulating bone loss and for determining risk of fracture in trabecular bone. Due to the computational efficiency, methods and products of the disclosed subject matter can be run on any suitable computing device, such as a personal computer 100 shown in FIG. 28, a workstation, a supercomputer, and the like.

[0169] In certain embodiments, systems for simulating bone loss can include a memory 102 that stores a three-dimensional image of a trabecular bone; and a processor 104 that: determines a number and location of bone resorption sites; carries out simulated bone resorption; identifies resorption sites associated with permanent bone loss; carries out a simulated bone formation process near resorption sites that are not identified as resorption cavities associated with permanent bone loss; carries out a simulated secondary mineralization process wherein the density of the individual trabeculae of the trabecular bone are increased; carries out a simulated positive bone balance process wherein the additional bones are formed in the trabeculae of the trabecular bone; carries out a finite element analysis based on the simulated bone formation to determine one or more mechanical properties of at least one trabecula of the trabecular bone; carries out a finite element analysis based on the simulated secondary mineralization process and/or the simulated positive bone balance process to determine one or more mechanical properties of at least one trabecular of the trabecular bone; and evaluates risk of fracture of the trabecular bone based on the finite element analysis.

[0170] In other embodiments, systems for determining risk of fracture in trabecular bone can include a memory that stores a three-dimensional image of a trabecular bone; and a processor that: removes voxels from trabeculae portion of the image to obtain a skeletonized structure of the trabecular bone, where the skeletonized structure of the trabecular bone is represented by voxels classified according to at least one topological classification; adds reconstruction voxels to at least voxels of a first topological classification; classifies the reconstruction voxels according to at least one topological classification that is related to the first topological classification; repeats said adding of reconstruction voxels and said classifying of the reconstruction voxels until a reconstructed image substantially representing the three-dimensional image of the trabecular bone is obtained; carries out morphological analysis and/or finite element analysis on the reconstructed image substantially representing the three-dimensional image of the trabecular bone; and evaluates risk of fracture of the trabecular bone based on the morphological analysis and/or the finite element analysis. [0171] As shown in FIG. 28, systems of the disclosed subject matter can further include a display device 106, an input device 108, an interface 110 to connect to an image capturing device, such as a μCT or a μMRI, and any other suitable devices (not shown). For example, systems of the disclosed subject matter can include a network interface to connect to the Internet, wherein the results of the analysis can be shared among various different users.

[0172] All patents, patent applications, and publications cited herein are hereby incorporated by reference herein in their entirety. The disclosures of these publications in their references are hereby incorporated by reference into this application in order to more fully describe the state of the art as known to those skilled therein as of the date of the invention described and claimed therein.

[0173] Upon review of the description and embodiments of the disclosed subject matter, those skilled in the art will understand that modifications and equivalent substitutions can be performed in carrying out the invention without departing from the essence of the invention. Thus, the invention is not meant to be limited by the embodiments and examples described explicitly above, but rather only by the claims which follow.

Claims

CLAIMSWhat is claimed is:

1. A method for simulating bone loss, the method comprising: reading a three-dimensional image of a trabecular bone; determining a number and location of bone resorption sites; and carrying out simulated bone resorption.

2. The method of claim 1, further comprising: identifying resorption sites associated with permanent bone loss; and carrying out a simulated bone formation process near resorption sites that are not identified as resorption cavities associated with permanent bone loss.

3. The method of claim 2, wherein the finite element analysis determines at least stress of the at least one trabecula; and the simulated bone formation is carried out near resorption sites that are in proximity to trabeculae with a stress level that is greater than a predetermined stress level.

4. The method of claim 2, further comprising: carrying out a finite element analysis based on the simulated bone formation to determine one or more mechanical properties of at least one trabecula of the trabecular bone.

5. The method of claim 4, further comprising: evaluating risk of fracture of the trabecular bone based on the finite element analysis.

6. The method of claim 2, further comprising at least one of: carrying out a simulated secondary mineralization process wherein the density of the individual trabeculae of the trabecular bone are increased; and carrying out a simulated positive bone balance process wherein the additional bones are formed in the trabeculae of the trabecular bone.

7. The method of claim 6, further comprising: carrying out a finite element analysis based on the simulated secondary mineralization process and/or the simulated positive bone balance process to determine one or more mechanical properties of at least one trabecular of the trabecular bone.

8. The method of claim 7, further comprising: evaluating risk of fracture of the trabecular bone based on the finite element analysis.

9. The method of claim 1, wherein the three-dimensional image of the trabecular bone is obtained from a micro-computed tomography scan or a micro magnetic resonance imaging scan.

10. A computer program product for simulating bone loss residing on a computer readable medium, the computer program product comprising instructions for causing a computer to: read a three-dimensional image of a trabecular bone; determine a number and location of bone resorption sites; and carry out simulated bone resorption.

11. The computer program product of claim 10, further comprising instructions for causing a computer to: identify resorption sites associated with permanent bone loss; carry out a simulated bone formation process near resorption sites that are not identified as resorption cavities associated with permanent bone loss; carry out a finite element analysis based on the simulated bone formation process to determine one or more mechanical properties of at least one trabecula of the trabecular bone; evaluate risk of fracture of the trabecular bone based on the finite element analysis.

12. A method for determining risk of fracture in trabecular bone, the method comprising: reading a three-dimensional image of a trabecular bone; removing voxels from a trabeculae portion of the image to obtain a skeletonized structure of the trabecular bone, where the skeletonized structure of the trabecular bone is represented by voxels classified according to at least one topological classification; and adding reconstruction voxels to at least voxels of a first topological classification.

13. The method of claim 12, further comprising: classifying the reconstruction voxels according to at least one topological classification that is related to the first topological classification.

14. The method of claim 13, further comprising: repeating said adding reconstruction voxels and said classifying the reconstruction voxels until a reconstructed image substantially representing the three-dimensional image of the trabecular bone is obtained.

15. The method of claim 14, further comprising: carrying out morphological analysis and/or finite element analysis on the reconstructed image substantially representing the three-dimensional image of the trabecular bone.

16. The method of claim 15, further comprising: evaluating risk of fracture of the trabecular bone based on the morphological analysis and/or the finite element analysis.

17. The method of claim 12, wherein the three-dimensional image of the trabecular bone is obtained from a micro -computed tomography scan or a micro magnetic resonance imaging scan.

18. A computer program product for determining risk of fracture in trabecular bone residing on a computer readable medium, the computer program product comprising instructions for causing a computer to: read a three-dimensional image of a trabecular bone; remove voxels from trabeculae portion of the image to obtain a skeletonized structure of the trabecular bone, where the skeletonized structure of the trabecular bone is represented by voxels classified according to at least one topological classification; and add reconstruction voxels to at least voxels of a first topological classification.

19. The computer program product of claim 18, further comprising instructions for causing a computer to: classify the reconstruction voxels according to at least one topological classifications that is related to the first topological classification; repeat said instruction to add reconstruction voxels and said instructions to classify the reconstruction voxels until a reconstructed image substantially representing the three- dimensional image of the trabecular bone is obtained; carry out morphological analysis and/or finite element analysis on the reconstructed image substantially representing the three-dimensional image of the trabecular bone; evaluate risk of fracture of the trabecular bone based on the morphological analysis and/or the finite element analysis.

20. A system for simulating bone loss, the system comprising: means for reading a three-dimensional image of a trabecular bone; means for determining a number and location of bone resorption sites; means for carrying out simulated bone resorption; means for identifying resorption sites associated with permanent bone loss; means for carrying out a simulated bone formation process near resorption sites that are not identified as resorption cavities associated with permanent bone loss; means for carrying out a simulated secondary mineralization process wherein the density of the individual trabeculae of the trabecular bone are increased; means for carrying out a simulated positive bone balance process wherein the additional bones are formed in the trabeculae of the trabecular bone; means for carrying out a finite element analysis based on the simulated bone formation to determine one or more mechanical properties of at least one trabecula of the trabecular bone; means for carrying out a finite element analysis based on the simulated secondary mineralization process and/or the simulated positive bone balance process to determine one or more mechanical properties of at least one trabecular of the trabecular bone; and means for evaluating risk of fracture of the trabecular bone based on the finite element analysis.

21. A system for determining risk of fracture in trabecular bone, the system comprising: means for reading a three-dimensional image of a trabecular bone; means for removing voxels from trabeculae portion of the image to obtain a skeletonized structure of the trabecular bone, where the skeletonized structure of the trabecular bone is represented by voxels classified according to at least one topological classification; means for adding reconstruction voxels to at least voxels of a first topological classification; means for classifying the reconstruction voxels according to at least one topological classification that is related to the first topological classification; means for repeating said adding reconstruction voxels and said classifying the reconstruction voxels until a reconstructed image substantially representing the three- dimensional image of the trabecular bone is obtained; means for carrying out morphological analysis and/or finite element analysis on the reconstructed image substantially representing the three-dimensional image of the trabecular bone; and means for evaluating risk of fracture of the trabecular bone based- on the morphological analysis and/or the finite element analysis.

22. A system for simulating bone loss, the system comprising: memory that stores a three-dimensional image of a trabecular bone; and a processor that: determines a number and location of bone resorption sites; carries out simulated bone resorption; identifies resorption sites associated with permanent bone loss; carries out a simulated bone formation process near resorption sites that are not identified as resorption cavities associated with permanent bone loss; carries out a simulated secondary mineralization process wherein the density of the individual trabeculae of the trabecular bone are increased; carries out a simulated positive bone balance process wherein the additional bones are formed in the trabeculae of the trabecular bone; carries out a finite element analysis based on the simulated bone formation to determine one or more mechanical properties of at least one trabecula of the trabecular bone; carries out a finite element analysis based on the simulated secondary mineralization process and/or the simulated positive bone balance process to determine one or more mechanical properties of at least one trabecular of the trabecular bone; and evaluates risk of fracture of the trabecular bone based on the finite element analysis.

23. A system for determining risk of fracture in trabecular bone, the system comprising: memory that stores a three-dimensional image of a trabecular bone; and a processor that: removes voxels from trabeculae portion of the image to obtain a skeletonized structure of the trabecular bone, where the skeletonized structure of the trabecular bone is represented by voxels classified according to at least one topological classification; adds reconstruction voxels to at least voxels of a first topological classification; classifies the reconstruction voxels according to at least one topological classification that is related to the first topological classification; repeats said adding of reconstruction voxels and said classifying of the reconstruction voxels until a reconstructed image substantially representing the three- dimensional image of the trabecular bone is obtained; carries out morphological analysis and/or finite element analysis on the reconstructed image substantially representing the three-dimensional image of the trabecular bone; and evaluates risk of fracture of the trabecular bone based on the morphological analysis and/or the finite element analysis.