- Email: [email protected]
Numerical simulation of the dependence of quantitative ultrasonic parameters on trabecular bone microarchitecture and elastic constants
Ultrasonics 44 (2006) e289–e294 www.elsevier.com/locate/ultras
Numerical simulation of the dependence of quantitative ultrasonic parameters on trabecular bone microarchitecture and elastic constants G. Haı¨at a, F. Padilla b
a,*
, R. Barkmann b, C.-C. Gluer b, P. Laugier
a
a Laboratoire d’Imagerie Parame´trique, CNRS UMR 7623 – Universite´ Paris 6, 15 rue de l’Ecole de Medecine, 75006 Paris, France Medizinische Physik, Klinik fu¨r Diagnostische Radiologie, Universita¨tsklinikum Schleswig-Holstein, Michaelisstraße 9, D-24105 Kiel, Germany
Available online 30 June 2006
Abstract Finite-difference numerical simulation of ultrasound propagation in complex media such as cancellous bone represents a fertile alternative to analytical approaches because it can manage the complex 3D bone structure by coupling the numerical computation with 3D numerical models of bone microarchitecture obtained from high-resolution imaging modalities. The objective of this work was to assess in silico the sensitivity of ultrasound parameters to controlled changes of microarchitecture and variation of elastic constants. The simulation software uses a finite-difference approach based on the Virieux numerical scheme. An incident plane wave was propagated through a volume of bone of approximately 5 · 5 · 8 mm3. The volumes were reconstructed from high-resolution micro-computed tomography data. An iterative numerical scenario of ‘‘virtual osteoporosis’’ was implemented using a dedicated image processing algorithm in order to modify the initial 3D microstructures. Numerical computations of wave propagation were performed at each step of the process. The sensitivity to bone material properties was also tested by changing the elastic constants of bone tissue. Our results suggest that ultrasonic variables (slope of the frequency-dependent attenuation coefficient and speed of sound) are mostly influenced by bone volume fraction. However, material properties and structure also appear to play a role. The impact of modifications of the stiffness coefficients remained lower than the variability caused by structural variations. This study emphasizes the potential of numerical computations tools coupled to realistic 3D structures to elucidate the physical mechanisms of interaction between ultrasound and bone structure and to assess the sensitivity of ultrasound variables to different bone properties. 2006 Elsevier B.V. All rights reserved. Keywords: Quantitative ultrasound; Bone; Virtual osteoporosis; Finite-difference; Porous media
1. Introduction In recent years, non-invasive quantitative ultrasound (QUS) techniques have become important modalities for the investigation of bone strength [1–5]. The measurement of ultrasonic properties of bone such as the speed of sound (SOS) or the slope of frequency-dependent attenuation (broadband ultrasonic attenuation BUA) using the transmission of elastic wave through the skeleton has been
*
Corresponding author. Tel.: +33 1 44 41 49 72; fax: +33 1 46 33 56 73. E-mail address: [email protected] (F. Padilla).
0041-624X/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ultras.2006.06.015
shown to convey important information related to bone strength [6,7]. In some studies, QUS parameters have been found to be closely related to site-matched estimates of bone density [8–11], a strong predictor of bone strength. Other reports point the ability of QUS parameters to reflect other non density-related bone characteristics [12–15]. QUS measurements of cancellous bone (i.e., the calcaneus) now represent a well established clinical modality for the prediction of fracture risk [16]. Nevertheless, despite their routine clinical use, the physics of ultrasound propagation through cancellous bone still remains unclear and the physical meaning of QUS parameters is still debated. This has led many clinicians to mistrust the technique.
e290
G. Haı¨at et al. / Ultrasonics 44 (2006) e289–e294
In the past, numerous aspects of ultrasonic bone characterization have been investigated and several theoretical models have been applied to the propagation of ultrasound in the MHz range through cancellous bone, including Biot [17–19], Schoenberg [20–22] and scattering models [19–23]. However, none of these approaches could account for all reported experimental results. In particular, the sensitivity of QUS measurement to several important microarchitectural characteristics of bone remains to be determined. Progress in this field has been hampered by the extreme complexity of the microarchitecture of bone. Trabecular bone consists of an anisotropic, heterogeneous and two phase medium made of a network of interconnected elastic rods or plates (the trabeculae) filled with a viscous fluid-like medium (the marrow). The difficulty of modelling the interaction of ultrasound propagation with trabecular bone arises from the fact that the typical size of the structural units (trabeculae, pores) is of the same order than the ultrasonic wave length at clinical frequencies (around 1.5 mm at 1 MHz). Developing a theory analytically including scattering, mode conversion and visco-elastic coupling rapidly become intractable. Finite-difference numerical simulation represents a promising alternative to analytical approaches because it can manage the complex three-dimensional (3D) bone structure by coupling the numerical computation with 3D numerical models of bone microarchitecture obtained from high-resolution imaging modalities [24]. Unlike former analytical approaches (Biot, Schonberg or scattering models) this numerical approach has been shown to be able to predict the major results observed in experiments [25], including the linear variation of the attenuation coefficient with frequency, negative velocity dispersion, existence of two compression waves (fast and slow waves), strong correlation of BUA and SOS to bone volume fraction. The objective of this work was to assess in silico the sensitivity of QUS parameters to controlled changes of microarchitecture and variation of elastic constants.
a standard desktop PC while maintaining a reasonably accurate representation of the trabecular structure. The bone volume fraction was defined as the total bone volume normalized to the total specimen volume (BV/TV). The microarchitecture was modified using a dedicated algorithm. Two scenarios of microarchitecture alteration were considered. The first one corresponds to an erosion of the trabecular network (i.e., decrease of BV/TV) and is assumed to model the thinning of the trabeculae which has been described to occur in to osteoporosis [27]. The second one corresponds to the complementary scenario in which a dilatation procedure is applied to the network (augmentation of BV/TV). The first step of the erosion method consisted in determining which voxels of the original trabecular network were adjacent to water. Fig. 1(a) shows a 2D image of the trabecular network on the XY plane perpendicular to the direction of ultrasound propagation. In the second step a randomly distributed erosion procedure was applied to those voxels. For each voxel, a random number n was chosen between 0 and 100. The considered voxel was suppressed only if n is inferior to x (the probability of suppression). This random erosion was selected in order to obtain progressive variations of BV/ TV. This procedure was applied with a probability of suppression of 50% and then of 100% (Fig. 1(b)). It was then applied again to this new structure, as long as the BV/TV remains greater to 5% (typical lower limit for BV/TV values). A dilatation procedure is then similarly applied to the original network leading to an increase of BV/TV (Fig. 1(c)), with the constraint that the BV/TV remained lower than 25%. For each sample, this process resulted in a set of 5–6 new 3D trabecular structures corresponding to the original network. In order to investigate the effect of the structure on QUS parameters independently of BV/TV, a similar method was applied to each sample until 3D structures with a constant BV/TV of 13% was reached. Therefore, the probability of
2. Material and method 2.1. 3D numerical structure models The numerical computations were coupled to 3D numerical models of trabecular structures. Toward this end, cylindrical cores were machined in 32 trabecular human specimens (typically 10 mm in height and 8 mm in diameter). The 3D microarchitecture of these samples was imaged using synchrotron radiation micro-computed tomography SR-lCT [26] at the ESRF (European Synchrotron Radiation Facility, France), with a voxel size of 10 lm in the three directions of space. Parallelepipedic domains with dimensions 5.6 · 5.6 mm2 in the transverse dimensions (X,Y) and 10 mm along the Z axis were extracted from each synchrotron acquisitions. From the initial resolution of 10 lm, the volumes were downsampled to a voxel size of 30 lm to match memory requirements on
Fig. 1. Illustration of the microarchitecture alteration scenarios. Images represent the intersection of the 3D trabecular structures with a plane perpendicular to the direction of propagation. Voxels corresponding to bone and to water are white and black saturated respectively. (a): Original sample, (b): eroded structure (probability of suppression of 100%), and (c): dilated cartographies (probability of addition of 100%).
G. Haı¨at et al. / Ultrasonics 44 (2006) e289–e294
voxel suppression or addition was adjusted in order to obtain structures with the same porosity for all 32 samples. 2.2. Simulation of the ultrasonic propagation Simulations of wave propagation through trabecular bone were performed using a software developed by our group (SimSonic), based on a finite-difference time-domain (FDTD) algorithm [28]. A numerical solution to the 3D linear elastic wave propagation equation is computed in both fluid and solid materials. Perfectly matched layers [29] have been implemented, avoiding unphysical reflections from the boundaries of the simulation mesh. The geometry of the propagation media is defined by a simple 3D regular mapping of the mechanical properties of the two media (water and bone). Complex phenomena such as reflection, refraction and mode conversion at the numerous bone interfaces are fully taken into account by the simulation tool. However, it does not account for absorption phenomena such as viscous effects. To simulate the propagation of plane waves through samples of limited transverse dimensions, symmetric boundary conditions were applied on the ZX and YZ faces of the simulation domain. The ultrasonic pressure source was defined over the whole XY plane at z = 0 mm. Therefore, propagation of plane waves through trabecular structures periodized by the boundary conditions in the X and Y directions were modelled by the simulation. A 5.6 · 5.6 mm2 plane receiver was set at z = 10 mm, recording the pressure averaged over its whole surface. The source emits a broadband pressure pulse centered at 1 MHz. For all computations, the interstitial phase was considered to be non viscous water, with density of 1 g/cm3 and acoustic wave velocity of 1500 m/s. Bone tissue was assumed to be isotropic and non absorbing, with a constant density of 1.85 g/cm3. To test the sensitivity of ultrasound parameters to bone material properties, five sets of stiffness coefficients indicated in Table 1 have been successively used to model bone. The reference properties correspond to a longitudinal velocity (VL) of 4000 m/s and transverse wave velocity (VT) of 1800 m/s, which are the typical elastic characteristics of cortical bone found in the literature [30]. The values of VL and VT were modified to correspond to a variation of ±10% which is assumed to stay within the physiological range.
e291
SOS and BUA were computed for each simulation. To derive the frequency-dependent attenuation, the ratio of the spectrum of the signal transmitted through bone to that of the reference signal was computed, the attenuation coefficient was deduced from the magnitude of the spectra ratio [31]. Here, SOS was computed with the first zero crossing method [32]. 3. Results Fig. 2 displays the variation of BUA (Fig. 2(a)) and SOS (Fig. 2(b)) as a function the BV/TV of the altered structures obtained from the erosion or dilatation procedure. Both parameters exhibit a strong positive correlation with BVTV. BUA and SOS were found to increase at an approximate rate of 2 dB/MHz/cm and 8 m/s respectively per % increase of BV/TV. Note that the dispersion of QUS variables for a given BV/TV value seems to be an increasing function of BV/TV. Mean value, standard deviation (SD) and range of variation of BUA and SOS obtained from simulations performed with the 32
Table 1 Elastic constants and corresponding wave propagation velocities of bone for the five set of parameters used in the simulation of ultrasonic propagation (density = 1.85 g/cm2)
Reference 1 2 3 4
C11 (GPa)
C44 (GPa)
VL (m/s)
VS (m/s)
29.6 23.976 35.816 29.6 29.6
5.994 5.994 5.994 4.736 7.4
4000 3600 4400 4000 4000
1800 1800 1800 1600 2000
Fig. 2. Variations of BUA (a) and SOS (b) as a function of BV/TV obtained from numerical simulations of ultrasonic propagation through in silico altered 3D structures. QUS parameters from the original (unaltered) structures are indicated by black dots and from the alteration process by crosses.
e292
G. Haı¨at et al. / Ultrasonics 44 (2006) e289–e294
Table 2 Mean, standard deviation (SD) and range of variation of SOS and BUA values obtained from simulations performed with 32 structures having the same BV/TV value (13%)
Mean ± SD Minimum Maximum
SOS (m/s)
BUA (dB/MHz/cm)
1586 ± 9 1567 1603
16.8 ± 1.95 13.32 21.3
All structures were derived from the 32 original structures.
structures of similar BV/TV (13%) are summarized in Table 2. These data indicate the variability (given by 4 SD) of QUS variables, which is due specifically to structure independently from bone quantity and material properties. For all 32 samples, the relative variations of BUA with the longitudinal and transverse velocities of bone, expressed in percentage, were computed as: BUAð2Þ BUAð1Þ BUAðrefÞ BUAð4Þ BUAð3Þ ¼ 100: BUAðrefÞ
where BUA(i) and BUA(Ref) are BUA obtained with the set of parameters #i and with the reference set of parameters. Similarly, the relative variations of SOS with VL and VT (noted DSOS_VL and DSOS_VT) were computed. These quantities represent the relative variation of QUS parameters for a relative variation of ±10% in the ultrasonic velocities of bone. DBUA_VL was consistently found to be positive, whereas for a few specimens (N = 4), negative values of DBUA_VT were found. In average, an increase of 20% in VL and in VT results in a relative increase of BUA of nearly 10%, with a relatively low standard deviation of ±2% for the change in VL and a much larger dispersion for VT (standard deviation of ±6.8%). The effects of VL and of VT were found to be independent of BV/TV (Fig. 3(a)). Similarly, DSOS_VL and DSOS_VT were consistently found to be positive (except DSOS_VT for one sample). In contrast to BUA, DSOS_VL and DSOS_VT were strongly related to BV/TV with a quasi linear variation (Fig. 3(b)).
DBUAV L ¼ 100: DBUAV L
ð1Þ
Fig. 3. Relative variations of BUA (a) and SOS (b) defined in Eq. (1) as a function of BV/TV. The relative variations correspond to variations in longitudinal (circles) and transverse (crosses) ultrasonic velocities of the bone matrix.
4. Discussion To our knowledge, this is the first report on the effect of virtual changes of bone quality (structural and material properties) on QUS parameters. We observed that changes in the trabecular structure has a direct impact on BUA. In this study, the simulation model does not account for absorption. and attenuation is assumed to be only due to scattering. Attenuation is thus a function of the size and shape of the scatterers, of the frequency as well as of the impedance mismatch between bone and water [33]. The scenario applied for the modification of the structure results in a thinning or thickening of the trabeculae. As we assumed that attenuation is due to scattering, we expected to observe a decrease of attenuation when the trabeculae are thinned, because their scattering crosssection is reduced, and less energy is thus scattered by the trabecular structure. The modifications of material properties also modify the scattering cross-section of the scatterers. As for the effect of scatterers size, an increase of the impedance mismatch between water and solid material is expected to produce an increase in the scattering cross-section, and consequently an increase of attenuation, as observed in the simulation. The dependence of BUA on VT is comparable in magnitude to the dependence of BUA on VL. It suggests that mode conversions and transverse wave propagation play an important role in ultrasonic attenuation in trabecular bone. The former arguments could only be qualitative, because as mentioned earlier, the exact mechanisms implied in the interaction between bone and ultrasound still remain unclear. In particular, further work is needed to understand why changes in VT induce a greater variability of BUA than changes in VL, or why the frequency dependence of the attenuation remains linear in specimens with different bone volume fraction and/or material properties.
G. Haı¨at et al. / Ultrasonics 44 (2006) e289–e294
e293
A simple homogenization model was developed in order to understand the behavior of SOS variations with BV/TV and material properties. Trabecular bone is modeled by a slice of homogeneous bone (in which only longitudinal waves propagates) and a slice of water, their respective dimension being determined by BV/TV. Hence, SOS writes as:
frequency-dependent attenuation coefficient and speed of sound) are mostly influenced by bone volume fraction. However, material properties and structure also appear to play a role. The impact of modifications of the stiffness coefficients remained lower than the variability caused by structural variations.
1 BV=TV ð1 BV=TVÞ ¼ þ SOS VL Vw
Acknowledgement ð2Þ
where Vw is the velocity in water. Changing the trabeculae thickness by erosion (resp. dilatation) results in a decrease (resp. increase) in SOS (see Fig. 2(b)). This general trend is found in the homogenization model (plotted with a dashed line), although it overestimates SOS values. This overestimation may be due to the presence of scattering (increase of wave path length) and/or to the influence of the transverse wave mode propagation. Moreover, the quasi-linear dependence of DSOS_VL with BV/TV (Fig. 3(b)) is in qualitative good agreement with our model. For low porosity values, this model overestimates the results obtained from simulation but a good agreement can be noticed for high porosity values. SOS was also found to be an increasing function of VT for all samples except for four. A significant (but weaker) linear dependence can also be observed between DSOS_VT and BV/TV. This aspect cannot be explained by the previous model and more work is needed to understand this behavior. As for BUA, the dependence of SOS on VT suggests that mode conversions and transverse wave propagation play an important role in ultrasonic propagation through trabecular bone. Finally, this study shows that the main determinant of QUS parameters appears to be BV/TV, as a variations of 20% of BV/TV induces a variation of about 45 dB/MHz/ cm for BUA and of about 250 m/s for SOS. For both BUA and SOS, the relative changes due to modifications of the stiffness coefficients (1.7 dB/MHz/cm for BUA and 15 m/s for SOS) remain much lower than the variability caused by structure variation for a given BV/TV of 13% (7.8 dB/MHz/cm for BUA and 36 m/s for SOS, figures given by 4 times the standard deviation indicated Table 2). Hence, the effect on QUS variables of microstructure variability appears to be more important than the effect of material properties. This variability may be due to different distributions of size and/or spacing of the scatterers between different structures. However, further studies are needed in order to clarify this last point. 5. Conclusion This study emphasizes the potential of numerical computations tools coupled to realistic 3D structures to elucidate the physical mechanisms of interaction between ultrasound and bone structure and to assess the sensitivity of ultrasound variables to different bone properties. Our results suggest that ultrasonic variables (slope of the
This work was supported by the European Commission (Contract no QLK6-CT-2002-02710). The authors would like to acknowledge Emmanuel Bossy for his large contribution to the numerical computations. References [1] E. Bossy, M. Talmant, P. Laugier, Bi-directional axial transmission can improve accuracy and precision of ultrasonic velocity measurement in cortical bone: a validation on test material, IEEE T. Ultrason., Ferr. 51 (1) (2004) 71–79. [2] P. Laugier et al., Ultrasound attenuation imaging in the os calcis: An improved method, Ultrason. Imag. 16 (2) (1994) 65–76. [3] C.M. Langton, S.B. Palmer, S.W. Porter, The measurement of broadband ultrasonic attenuation in cancellous bone, Eng. Med. 13 (2) (1984) 89–91. [4] R. Barkmann et al., A new method for quantitative ultrasound measurements at multiple skeletal sites, J. Clin. Densit. 3 (2000) 1–7. [5] R. Barkmann et al., Assessment of the geometry of human finger phalanges using quantitative ultrasound in vivo, Osteoporos Int. 11 (9) (2000) 745–755. [6] C.F. Njeh et al., Prediction of human femoral bone strength using ultrasound velocity and BMD: an in vitro study, Osteoporosis Int. 7 (1997) 471–477. [7] M.L. Bouxsein, S.E. Radloff, Quantitative ultrasound of the calcaneus reflects the mechanical properties of calcaneal trabecular bone, J. Bone Miner. Rea. 12 (5) (1997) 839–846. [8] C. Chappard et al., Assessment of the relationship between broadband ultrasound attenuation and bone mineral density at the calcaneus using BUA imaging and DXA, Osteoporosis Int. 7 (4) (1997) 316–322. [9] F. Padilla et al., In vitro ultrasound measurement at the human femur, Calcif. Tissue. Int. 75 (5) (2004) 421–430. [10] G. Haı¨at et al., Optimal prediction of bone mineral density with ultrasonic measurements in excised human femur, Calcif. Tissue. Int. 77 (3) (2005) 186–192. [11] G. Haı¨at et al., In vitro speed of sound measurement at intact human femur specimens, Ultrasound. Med. Biol. 31 (7) (2005) 987–996. [12] C.C. Gluer, C.Y. Wu, H.K. Genant, Broadband ultrasound attenuation signals depends on trabecular orientation: an in vitro study, Osteoporosis Int. 3 (1993) 185–191. [13] P.H.F. Nicholson, M.J. Haddaway, M.W.J. Davie, The dependence of ultrasonic properties on orientation in human vertebral bone, Phys. Med. Biol. 39 (1994) 1013–1024. [14] D. Hans et al., Ultrasound velocity of trabecular cubes reflects mainly bone density and elasticity, Calcif. Tissue. Int. 64 (1999) 18– 23. [15] K.A. Wear, Anisotropy of ultrasonic backscatter and attenuation from human calcaneus : Implications for relative roles of absorption and scattering in determining attenuation, J. Acoust. Soc. Am. 107 (6) (2000) 3474–3479. [16] D. Hans et al., Ultrasonographic heel measurements to predict hip fracture in elderly women: the EPIDOS prospective study, Lancet 348 (9026) (1996) 511–514. [17] A. Hosokawa, T. Otani, Ultrasonic wave propagation in bovine cancellous bone, J. Acoust. Soc. Am. 101 (1) (1997) 1–5.
e294
G. Haı¨at et al. / Ultrasonics 44 (2006) e289–e294
[18] Z.E. Fellah et al., Ultrasonic wave propagation in human cancellous bone: Application of Biot theory, J. Acoust. Soc. Am. 116 (1) (2004) 61–73. [19] F. Luppe, J.M. Conoir, H. Franklin, Scattering by a fluid cylinder in a porous medium: Application to trabecular bone, J. Acoust. Soc. Am. 111 (6) (2002) 2573–2582. [20] E.R. Hughes et al., Ultrasonic propagation in cancellous bone : A new stratified model, Ultrasound Med. Biol. 25 (5) (1999) 811–821. [21] F. Padilla, P. Laugier, Phase and group velocities of fast and slow compressional waves in trabecular bone, J. Acoust. Soc. Am. 108 (4) (2000) 1949–1952. [22] K.A. Wear, A stratified model to predict dispersion in trabecular bone, IEEE T. Ultrason. Ferr. 48 (4) (2001) 1079–1083. [23] P.H.F. Nicholson et al., Scattering of ultrasound in cancellous bone:predictions from a theoretical model, J. Biomech. 33 (2000) 503–506. [24] G. Luo et al., Computational methods for ultrasonic bone assessment, Ultrasound Med. Biol. 25 (5) (1999) 823–830. [25] E. Bossy, F. Padilla, P. Laugier, Three-dimensional simulation of ultrasound propagation through trabecular bone structures measured by synchrotron microtomography, Phys. Med. Biol. 50 (23) (2005) 5545–5556.
[26] M. Salome et al., A synchrotron radiation microtomography system for the analysis of trabecular bone samples, Med. Phys. 26 (10) (1999) 2194–2204. [27] L. Pothuaud et al., Fractal dimension of trabecular bone projection texture is related to three-dimensional microarchitecture, J. Bone Miner. Res. 15 (4) (2000) 691–699. [28] E. Bossy, M. Talmant, P. Laugier, Three-dimensional simulations of ultrasonic axial transmission velocity measurement on cortical bone models, J. Acoust. Soc. Am. 115 (5 Pt 1) (2004) 2314–2324. [29] F. Collino, C. Tsogka, Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media, Geophysics 66 (1) (2001) 294–307. [30] C.F. Njeh et al., Ultrasonic Studies of Cortical Bone In vitro, in Quantitative Ultrasound Assessment of Osteoporosis and Bone Status, Martin Dunitz, London, 1999, p. 177. [31] P. Droin, G. Berger, P. Laugier, Velocity dispersion of acoustic waves in cancellous bone, IEEE T. Ultrason. Ferr. 45 (3) (1998) 581–592. [32] G. Haı¨at et al., Effects of frequency-dependent attenuation and velocity dispersion on in vitro ultrasound velocity measurements in intact human femur specimens, IEEE T. Ultrason Ferr. 53 (1) (2006) 39–51. [33] P. Morse, K. Ingard, Theoretical Acoustics, Princetown University Press, Pricetown, NJ, 1986.
Numerical simulation of the dependence of quantitative ultrasonic parameters on trabecular bone microarchitecture and elastic constants G. Haı¨at a, F. Padilla b
a,*
, R. Barkmann b, C.-C. Gluer b, P. Laugier
a
a Laboratoire d’Imagerie Parame´trique, CNRS UMR 7623 – Universite´ Paris 6, 15 rue de l’Ecole de Medecine, 75006 Paris, France Medizinische Physik, Klinik fu¨r Diagnostische Radiologie, Universita¨tsklinikum Schleswig-Holstein, Michaelisstraße 9, D-24105 Kiel, Germany
Available online 30 June 2006
Abstract Finite-difference numerical simulation of ultrasound propagation in complex media such as cancellous bone represents a fertile alternative to analytical approaches because it can manage the complex 3D bone structure by coupling the numerical computation with 3D numerical models of bone microarchitecture obtained from high-resolution imaging modalities. The objective of this work was to assess in silico the sensitivity of ultrasound parameters to controlled changes of microarchitecture and variation of elastic constants. The simulation software uses a finite-difference approach based on the Virieux numerical scheme. An incident plane wave was propagated through a volume of bone of approximately 5 · 5 · 8 mm3. The volumes were reconstructed from high-resolution micro-computed tomography data. An iterative numerical scenario of ‘‘virtual osteoporosis’’ was implemented using a dedicated image processing algorithm in order to modify the initial 3D microstructures. Numerical computations of wave propagation were performed at each step of the process. The sensitivity to bone material properties was also tested by changing the elastic constants of bone tissue. Our results suggest that ultrasonic variables (slope of the frequency-dependent attenuation coefficient and speed of sound) are mostly influenced by bone volume fraction. However, material properties and structure also appear to play a role. The impact of modifications of the stiffness coefficients remained lower than the variability caused by structural variations. This study emphasizes the potential of numerical computations tools coupled to realistic 3D structures to elucidate the physical mechanisms of interaction between ultrasound and bone structure and to assess the sensitivity of ultrasound variables to different bone properties. 2006 Elsevier B.V. All rights reserved. Keywords: Quantitative ultrasound; Bone; Virtual osteoporosis; Finite-difference; Porous media
1. Introduction In recent years, non-invasive quantitative ultrasound (QUS) techniques have become important modalities for the investigation of bone strength [1–5]. The measurement of ultrasonic properties of bone such as the speed of sound (SOS) or the slope of frequency-dependent attenuation (broadband ultrasonic attenuation BUA) using the transmission of elastic wave through the skeleton has been
*
Corresponding author. Tel.: +33 1 44 41 49 72; fax: +33 1 46 33 56 73. E-mail address: [email protected] (F. Padilla).
0041-624X/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ultras.2006.06.015
shown to convey important information related to bone strength [6,7]. In some studies, QUS parameters have been found to be closely related to site-matched estimates of bone density [8–11], a strong predictor of bone strength. Other reports point the ability of QUS parameters to reflect other non density-related bone characteristics [12–15]. QUS measurements of cancellous bone (i.e., the calcaneus) now represent a well established clinical modality for the prediction of fracture risk [16]. Nevertheless, despite their routine clinical use, the physics of ultrasound propagation through cancellous bone still remains unclear and the physical meaning of QUS parameters is still debated. This has led many clinicians to mistrust the technique.
e290
G. Haı¨at et al. / Ultrasonics 44 (2006) e289–e294
In the past, numerous aspects of ultrasonic bone characterization have been investigated and several theoretical models have been applied to the propagation of ultrasound in the MHz range through cancellous bone, including Biot [17–19], Schoenberg [20–22] and scattering models [19–23]. However, none of these approaches could account for all reported experimental results. In particular, the sensitivity of QUS measurement to several important microarchitectural characteristics of bone remains to be determined. Progress in this field has been hampered by the extreme complexity of the microarchitecture of bone. Trabecular bone consists of an anisotropic, heterogeneous and two phase medium made of a network of interconnected elastic rods or plates (the trabeculae) filled with a viscous fluid-like medium (the marrow). The difficulty of modelling the interaction of ultrasound propagation with trabecular bone arises from the fact that the typical size of the structural units (trabeculae, pores) is of the same order than the ultrasonic wave length at clinical frequencies (around 1.5 mm at 1 MHz). Developing a theory analytically including scattering, mode conversion and visco-elastic coupling rapidly become intractable. Finite-difference numerical simulation represents a promising alternative to analytical approaches because it can manage the complex three-dimensional (3D) bone structure by coupling the numerical computation with 3D numerical models of bone microarchitecture obtained from high-resolution imaging modalities [24]. Unlike former analytical approaches (Biot, Schonberg or scattering models) this numerical approach has been shown to be able to predict the major results observed in experiments [25], including the linear variation of the attenuation coefficient with frequency, negative velocity dispersion, existence of two compression waves (fast and slow waves), strong correlation of BUA and SOS to bone volume fraction. The objective of this work was to assess in silico the sensitivity of QUS parameters to controlled changes of microarchitecture and variation of elastic constants.
a standard desktop PC while maintaining a reasonably accurate representation of the trabecular structure. The bone volume fraction was defined as the total bone volume normalized to the total specimen volume (BV/TV). The microarchitecture was modified using a dedicated algorithm. Two scenarios of microarchitecture alteration were considered. The first one corresponds to an erosion of the trabecular network (i.e., decrease of BV/TV) and is assumed to model the thinning of the trabeculae which has been described to occur in to osteoporosis [27]. The second one corresponds to the complementary scenario in which a dilatation procedure is applied to the network (augmentation of BV/TV). The first step of the erosion method consisted in determining which voxels of the original trabecular network were adjacent to water. Fig. 1(a) shows a 2D image of the trabecular network on the XY plane perpendicular to the direction of ultrasound propagation. In the second step a randomly distributed erosion procedure was applied to those voxels. For each voxel, a random number n was chosen between 0 and 100. The considered voxel was suppressed only if n is inferior to x (the probability of suppression). This random erosion was selected in order to obtain progressive variations of BV/ TV. This procedure was applied with a probability of suppression of 50% and then of 100% (Fig. 1(b)). It was then applied again to this new structure, as long as the BV/TV remains greater to 5% (typical lower limit for BV/TV values). A dilatation procedure is then similarly applied to the original network leading to an increase of BV/TV (Fig. 1(c)), with the constraint that the BV/TV remained lower than 25%. For each sample, this process resulted in a set of 5–6 new 3D trabecular structures corresponding to the original network. In order to investigate the effect of the structure on QUS parameters independently of BV/TV, a similar method was applied to each sample until 3D structures with a constant BV/TV of 13% was reached. Therefore, the probability of
2. Material and method 2.1. 3D numerical structure models The numerical computations were coupled to 3D numerical models of trabecular structures. Toward this end, cylindrical cores were machined in 32 trabecular human specimens (typically 10 mm in height and 8 mm in diameter). The 3D microarchitecture of these samples was imaged using synchrotron radiation micro-computed tomography SR-lCT [26] at the ESRF (European Synchrotron Radiation Facility, France), with a voxel size of 10 lm in the three directions of space. Parallelepipedic domains with dimensions 5.6 · 5.6 mm2 in the transverse dimensions (X,Y) and 10 mm along the Z axis were extracted from each synchrotron acquisitions. From the initial resolution of 10 lm, the volumes were downsampled to a voxel size of 30 lm to match memory requirements on
Fig. 1. Illustration of the microarchitecture alteration scenarios. Images represent the intersection of the 3D trabecular structures with a plane perpendicular to the direction of propagation. Voxels corresponding to bone and to water are white and black saturated respectively. (a): Original sample, (b): eroded structure (probability of suppression of 100%), and (c): dilated cartographies (probability of addition of 100%).
G. Haı¨at et al. / Ultrasonics 44 (2006) e289–e294
voxel suppression or addition was adjusted in order to obtain structures with the same porosity for all 32 samples. 2.2. Simulation of the ultrasonic propagation Simulations of wave propagation through trabecular bone were performed using a software developed by our group (SimSonic), based on a finite-difference time-domain (FDTD) algorithm [28]. A numerical solution to the 3D linear elastic wave propagation equation is computed in both fluid and solid materials. Perfectly matched layers [29] have been implemented, avoiding unphysical reflections from the boundaries of the simulation mesh. The geometry of the propagation media is defined by a simple 3D regular mapping of the mechanical properties of the two media (water and bone). Complex phenomena such as reflection, refraction and mode conversion at the numerous bone interfaces are fully taken into account by the simulation tool. However, it does not account for absorption phenomena such as viscous effects. To simulate the propagation of plane waves through samples of limited transverse dimensions, symmetric boundary conditions were applied on the ZX and YZ faces of the simulation domain. The ultrasonic pressure source was defined over the whole XY plane at z = 0 mm. Therefore, propagation of plane waves through trabecular structures periodized by the boundary conditions in the X and Y directions were modelled by the simulation. A 5.6 · 5.6 mm2 plane receiver was set at z = 10 mm, recording the pressure averaged over its whole surface. The source emits a broadband pressure pulse centered at 1 MHz. For all computations, the interstitial phase was considered to be non viscous water, with density of 1 g/cm3 and acoustic wave velocity of 1500 m/s. Bone tissue was assumed to be isotropic and non absorbing, with a constant density of 1.85 g/cm3. To test the sensitivity of ultrasound parameters to bone material properties, five sets of stiffness coefficients indicated in Table 1 have been successively used to model bone. The reference properties correspond to a longitudinal velocity (VL) of 4000 m/s and transverse wave velocity (VT) of 1800 m/s, which are the typical elastic characteristics of cortical bone found in the literature [30]. The values of VL and VT were modified to correspond to a variation of ±10% which is assumed to stay within the physiological range.
e291
SOS and BUA were computed for each simulation. To derive the frequency-dependent attenuation, the ratio of the spectrum of the signal transmitted through bone to that of the reference signal was computed, the attenuation coefficient was deduced from the magnitude of the spectra ratio [31]. Here, SOS was computed with the first zero crossing method [32]. 3. Results Fig. 2 displays the variation of BUA (Fig. 2(a)) and SOS (Fig. 2(b)) as a function the BV/TV of the altered structures obtained from the erosion or dilatation procedure. Both parameters exhibit a strong positive correlation with BVTV. BUA and SOS were found to increase at an approximate rate of 2 dB/MHz/cm and 8 m/s respectively per % increase of BV/TV. Note that the dispersion of QUS variables for a given BV/TV value seems to be an increasing function of BV/TV. Mean value, standard deviation (SD) and range of variation of BUA and SOS obtained from simulations performed with the 32
Table 1 Elastic constants and corresponding wave propagation velocities of bone for the five set of parameters used in the simulation of ultrasonic propagation (density = 1.85 g/cm2)
Reference 1 2 3 4
C11 (GPa)
C44 (GPa)
VL (m/s)
VS (m/s)
29.6 23.976 35.816 29.6 29.6
5.994 5.994 5.994 4.736 7.4
4000 3600 4400 4000 4000
1800 1800 1800 1600 2000
Fig. 2. Variations of BUA (a) and SOS (b) as a function of BV/TV obtained from numerical simulations of ultrasonic propagation through in silico altered 3D structures. QUS parameters from the original (unaltered) structures are indicated by black dots and from the alteration process by crosses.
e292
G. Haı¨at et al. / Ultrasonics 44 (2006) e289–e294
Table 2 Mean, standard deviation (SD) and range of variation of SOS and BUA values obtained from simulations performed with 32 structures having the same BV/TV value (13%)
Mean ± SD Minimum Maximum
SOS (m/s)
BUA (dB/MHz/cm)
1586 ± 9 1567 1603
16.8 ± 1.95 13.32 21.3
All structures were derived from the 32 original structures.
structures of similar BV/TV (13%) are summarized in Table 2. These data indicate the variability (given by 4 SD) of QUS variables, which is due specifically to structure independently from bone quantity and material properties. For all 32 samples, the relative variations of BUA with the longitudinal and transverse velocities of bone, expressed in percentage, were computed as: BUAð2Þ BUAð1Þ BUAðrefÞ BUAð4Þ BUAð3Þ ¼ 100: BUAðrefÞ
where BUA(i) and BUA(Ref) are BUA obtained with the set of parameters #i and with the reference set of parameters. Similarly, the relative variations of SOS with VL and VT (noted DSOS_VL and DSOS_VT) were computed. These quantities represent the relative variation of QUS parameters for a relative variation of ±10% in the ultrasonic velocities of bone. DBUA_VL was consistently found to be positive, whereas for a few specimens (N = 4), negative values of DBUA_VT were found. In average, an increase of 20% in VL and in VT results in a relative increase of BUA of nearly 10%, with a relatively low standard deviation of ±2% for the change in VL and a much larger dispersion for VT (standard deviation of ±6.8%). The effects of VL and of VT were found to be independent of BV/TV (Fig. 3(a)). Similarly, DSOS_VL and DSOS_VT were consistently found to be positive (except DSOS_VT for one sample). In contrast to BUA, DSOS_VL and DSOS_VT were strongly related to BV/TV with a quasi linear variation (Fig. 3(b)).
DBUAV L ¼ 100: DBUAV L
ð1Þ
Fig. 3. Relative variations of BUA (a) and SOS (b) defined in Eq. (1) as a function of BV/TV. The relative variations correspond to variations in longitudinal (circles) and transverse (crosses) ultrasonic velocities of the bone matrix.
4. Discussion To our knowledge, this is the first report on the effect of virtual changes of bone quality (structural and material properties) on QUS parameters. We observed that changes in the trabecular structure has a direct impact on BUA. In this study, the simulation model does not account for absorption. and attenuation is assumed to be only due to scattering. Attenuation is thus a function of the size and shape of the scatterers, of the frequency as well as of the impedance mismatch between bone and water [33]. The scenario applied for the modification of the structure results in a thinning or thickening of the trabeculae. As we assumed that attenuation is due to scattering, we expected to observe a decrease of attenuation when the trabeculae are thinned, because their scattering crosssection is reduced, and less energy is thus scattered by the trabecular structure. The modifications of material properties also modify the scattering cross-section of the scatterers. As for the effect of scatterers size, an increase of the impedance mismatch between water and solid material is expected to produce an increase in the scattering cross-section, and consequently an increase of attenuation, as observed in the simulation. The dependence of BUA on VT is comparable in magnitude to the dependence of BUA on VL. It suggests that mode conversions and transverse wave propagation play an important role in ultrasonic attenuation in trabecular bone. The former arguments could only be qualitative, because as mentioned earlier, the exact mechanisms implied in the interaction between bone and ultrasound still remain unclear. In particular, further work is needed to understand why changes in VT induce a greater variability of BUA than changes in VL, or why the frequency dependence of the attenuation remains linear in specimens with different bone volume fraction and/or material properties.
G. Haı¨at et al. / Ultrasonics 44 (2006) e289–e294
e293
A simple homogenization model was developed in order to understand the behavior of SOS variations with BV/TV and material properties. Trabecular bone is modeled by a slice of homogeneous bone (in which only longitudinal waves propagates) and a slice of water, their respective dimension being determined by BV/TV. Hence, SOS writes as:
frequency-dependent attenuation coefficient and speed of sound) are mostly influenced by bone volume fraction. However, material properties and structure also appear to play a role. The impact of modifications of the stiffness coefficients remained lower than the variability caused by structural variations.
1 BV=TV ð1 BV=TVÞ ¼ þ SOS VL Vw
Acknowledgement ð2Þ
where Vw is the velocity in water. Changing the trabeculae thickness by erosion (resp. dilatation) results in a decrease (resp. increase) in SOS (see Fig. 2(b)). This general trend is found in the homogenization model (plotted with a dashed line), although it overestimates SOS values. This overestimation may be due to the presence of scattering (increase of wave path length) and/or to the influence of the transverse wave mode propagation. Moreover, the quasi-linear dependence of DSOS_VL with BV/TV (Fig. 3(b)) is in qualitative good agreement with our model. For low porosity values, this model overestimates the results obtained from simulation but a good agreement can be noticed for high porosity values. SOS was also found to be an increasing function of VT for all samples except for four. A significant (but weaker) linear dependence can also be observed between DSOS_VT and BV/TV. This aspect cannot be explained by the previous model and more work is needed to understand this behavior. As for BUA, the dependence of SOS on VT suggests that mode conversions and transverse wave propagation play an important role in ultrasonic propagation through trabecular bone. Finally, this study shows that the main determinant of QUS parameters appears to be BV/TV, as a variations of 20% of BV/TV induces a variation of about 45 dB/MHz/ cm for BUA and of about 250 m/s for SOS. For both BUA and SOS, the relative changes due to modifications of the stiffness coefficients (1.7 dB/MHz/cm for BUA and 15 m/s for SOS) remain much lower than the variability caused by structure variation for a given BV/TV of 13% (7.8 dB/MHz/cm for BUA and 36 m/s for SOS, figures given by 4 times the standard deviation indicated Table 2). Hence, the effect on QUS variables of microstructure variability appears to be more important than the effect of material properties. This variability may be due to different distributions of size and/or spacing of the scatterers between different structures. However, further studies are needed in order to clarify this last point. 5. Conclusion This study emphasizes the potential of numerical computations tools coupled to realistic 3D structures to elucidate the physical mechanisms of interaction between ultrasound and bone structure and to assess the sensitivity of ultrasound variables to different bone properties. Our results suggest that ultrasonic variables (slope of the
This work was supported by the European Commission (Contract no QLK6-CT-2002-02710). The authors would like to acknowledge Emmanuel Bossy for his large contribution to the numerical computations. References [1] E. Bossy, M. Talmant, P. Laugier, Bi-directional axial transmission can improve accuracy and precision of ultrasonic velocity measurement in cortical bone: a validation on test material, IEEE T. Ultrason., Ferr. 51 (1) (2004) 71–79. [2] P. Laugier et al., Ultrasound attenuation imaging in the os calcis: An improved method, Ultrason. Imag. 16 (2) (1994) 65–76. [3] C.M. Langton, S.B. Palmer, S.W. Porter, The measurement of broadband ultrasonic attenuation in cancellous bone, Eng. Med. 13 (2) (1984) 89–91. [4] R. Barkmann et al., A new method for quantitative ultrasound measurements at multiple skeletal sites, J. Clin. Densit. 3 (2000) 1–7. [5] R. Barkmann et al., Assessment of the geometry of human finger phalanges using quantitative ultrasound in vivo, Osteoporos Int. 11 (9) (2000) 745–755. [6] C.F. Njeh et al., Prediction of human femoral bone strength using ultrasound velocity and BMD: an in vitro study, Osteoporosis Int. 7 (1997) 471–477. [7] M.L. Bouxsein, S.E. Radloff, Quantitative ultrasound of the calcaneus reflects the mechanical properties of calcaneal trabecular bone, J. Bone Miner. Rea. 12 (5) (1997) 839–846. [8] C. Chappard et al., Assessment of the relationship between broadband ultrasound attenuation and bone mineral density at the calcaneus using BUA imaging and DXA, Osteoporosis Int. 7 (4) (1997) 316–322. [9] F. Padilla et al., In vitro ultrasound measurement at the human femur, Calcif. Tissue. Int. 75 (5) (2004) 421–430. [10] G. Haı¨at et al., Optimal prediction of bone mineral density with ultrasonic measurements in excised human femur, Calcif. Tissue. Int. 77 (3) (2005) 186–192. [11] G. Haı¨at et al., In vitro speed of sound measurement at intact human femur specimens, Ultrasound. Med. Biol. 31 (7) (2005) 987–996. [12] C.C. Gluer, C.Y. Wu, H.K. Genant, Broadband ultrasound attenuation signals depends on trabecular orientation: an in vitro study, Osteoporosis Int. 3 (1993) 185–191. [13] P.H.F. Nicholson, M.J. Haddaway, M.W.J. Davie, The dependence of ultrasonic properties on orientation in human vertebral bone, Phys. Med. Biol. 39 (1994) 1013–1024. [14] D. Hans et al., Ultrasound velocity of trabecular cubes reflects mainly bone density and elasticity, Calcif. Tissue. Int. 64 (1999) 18– 23. [15] K.A. Wear, Anisotropy of ultrasonic backscatter and attenuation from human calcaneus : Implications for relative roles of absorption and scattering in determining attenuation, J. Acoust. Soc. Am. 107 (6) (2000) 3474–3479. [16] D. Hans et al., Ultrasonographic heel measurements to predict hip fracture in elderly women: the EPIDOS prospective study, Lancet 348 (9026) (1996) 511–514. [17] A. Hosokawa, T. Otani, Ultrasonic wave propagation in bovine cancellous bone, J. Acoust. Soc. Am. 101 (1) (1997) 1–5.
e294
G. Haı¨at et al. / Ultrasonics 44 (2006) e289–e294
[18] Z.E. Fellah et al., Ultrasonic wave propagation in human cancellous bone: Application of Biot theory, J. Acoust. Soc. Am. 116 (1) (2004) 61–73. [19] F. Luppe, J.M. Conoir, H. Franklin, Scattering by a fluid cylinder in a porous medium: Application to trabecular bone, J. Acoust. Soc. Am. 111 (6) (2002) 2573–2582. [20] E.R. Hughes et al., Ultrasonic propagation in cancellous bone : A new stratified model, Ultrasound Med. Biol. 25 (5) (1999) 811–821. [21] F. Padilla, P. Laugier, Phase and group velocities of fast and slow compressional waves in trabecular bone, J. Acoust. Soc. Am. 108 (4) (2000) 1949–1952. [22] K.A. Wear, A stratified model to predict dispersion in trabecular bone, IEEE T. Ultrason. Ferr. 48 (4) (2001) 1079–1083. [23] P.H.F. Nicholson et al., Scattering of ultrasound in cancellous bone:predictions from a theoretical model, J. Biomech. 33 (2000) 503–506. [24] G. Luo et al., Computational methods for ultrasonic bone assessment, Ultrasound Med. Biol. 25 (5) (1999) 823–830. [25] E. Bossy, F. Padilla, P. Laugier, Three-dimensional simulation of ultrasound propagation through trabecular bone structures measured by synchrotron microtomography, Phys. Med. Biol. 50 (23) (2005) 5545–5556.
[26] M. Salome et al., A synchrotron radiation microtomography system for the analysis of trabecular bone samples, Med. Phys. 26 (10) (1999) 2194–2204. [27] L. Pothuaud et al., Fractal dimension of trabecular bone projection texture is related to three-dimensional microarchitecture, J. Bone Miner. Res. 15 (4) (2000) 691–699. [28] E. Bossy, M. Talmant, P. Laugier, Three-dimensional simulations of ultrasonic axial transmission velocity measurement on cortical bone models, J. Acoust. Soc. Am. 115 (5 Pt 1) (2004) 2314–2324. [29] F. Collino, C. Tsogka, Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media, Geophysics 66 (1) (2001) 294–307. [30] C.F. Njeh et al., Ultrasonic Studies of Cortical Bone In vitro, in Quantitative Ultrasound Assessment of Osteoporosis and Bone Status, Martin Dunitz, London, 1999, p. 177. [31] P. Droin, G. Berger, P. Laugier, Velocity dispersion of acoustic waves in cancellous bone, IEEE T. Ultrason. Ferr. 45 (3) (1998) 581–592. [32] G. Haı¨at et al., Effects of frequency-dependent attenuation and velocity dispersion on in vitro ultrasound velocity measurements in intact human femur specimens, IEEE T. Ultrason Ferr. 53 (1) (2006) 39–51. [33] P. Morse, K. Ingard, Theoretical Acoustics, Princetown University Press, Pricetown, NJ, 1986.