- Email: [email protected]
Understanding velocity of sound in trabecular bone via computer simulations
Computers in Biology and Medicine 36 (2006) 439 – 447 www.intl.elsevierhealth.com/journals/cobm
Understanding velocity of sound in trabecular bone via computer simulations Raafat George Saadéa,∗ , George Tsoukasb , John Caminisc a Department of Decision Sciences and MIS, John Molson School of Business, Concordia University, 1455 De Maisonneuve
Building, Montreal, Queb., Canada H3G 1M8 b Center for Bone Metabolism, Faculty of Medicine, McGill University, Canada c Global Clinical Development & Medical Affairs, East Hanover, NJ 07936-1080, USA
Received 20 December 2004; received in revised form 28 February 2005; accepted 14 March 2005
Abstract Osteoporosis is a condition characterized by low bone mass and micro-architectural deterioration leading to non-traumatic fractures of the skeleton. It is a potentially debilitating condition especially for senior citizens. There is growing interest for quantitative ultrasound to measure bone mineral density. However, understanding of ultrasound–bone interaction is limited. Simulating ultrasound propagation of bone can help us better understand ultrasound–bone interaction and provide insight into bone architecture. In this study a mathematical model for the propagation of sound in bone is presented. Results demonstrate the suitability of the proposed modeling approach and the model’s capability to reproduce conditions in the lab. 䉷 2005 Published by Elsevier Ltd. Keywords: Artificial intelligence; Clinical; Disease management; Osteoporosis
1. Introduction Osteoporosis is a condition characterized by low bone mass and micro-architectural deterioration leading to non-traumatic fractures of the skeleton. It is a potentially debilitating condition especially for senior citizens and carries a direct cost of over 15 billion dollars a year in medical expenses, in the US alone [5]. Early diagnosis of osteoporosis is critical, if the risk of fracture is to be reduced by optimizing ∗ Corresponding author. Tel.: +1 450 682 5563; fax: +1 514 848 2424x2988.
E-mail address: [email protected] (R.G. Saadé). 0010-4825/$ - see front matter 䉷 2005 Published by Elsevier Ltd. doi:10.1016/j.compbiomed.2005.03.008
440
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
prevention and treatment strategies. Fracture is the principle outcome of osteoporosis, and it is essential to then to understand the determinants of fracture risk. One major determinant is intrinsic bone strength. Bone strength and its ability to resist fracture, is influenced by bone structure which includes bone size, shape, and micro-architectural characteristics [1]. Subsequently, and based on clinical observations; material properties of bone, namely density and elasticity, play an important role in risk fracture making them important surrogates for risk assessment. Bone from older individuals demonstrates decreased “tensile plastic deformation” (a result of compromised elastic properties of bone) that are attributed to molecular changes in either organic or the inorganic bone matrix [2]. This age-related reduction in the ability of bone to absorb and redistribute energy resulting from a fall may provide important information in assessing risk or diagnosing osteoporosis. Hence, the loss of energy absorption capability may be a significant factor increasing the risk of fracture in older women with low bone mass. Bone mineral density (BMD) is the most common and accepted osteoporosis diagnostic index. BMD is measured by dual energy X-ray absorptiometry, or DXA [3], and has been reported to explain about 60–75% of the variance in strength, and the remaining variance could be due to other factors such as accumulated fatigue damage, inferior bone micro-architecture and the state of bone remodeling [4]. However, the use of DXA has noted limitations: perception of excess patient exposed to radiation, expensive equipment, and its measurement of only apparent density, which is an image processing of different pixel densities, limiting its ability to proportionally correlate with strength [5]. Due to the inherent limitation of DXA, there is a growing interest in the use of quantitative ultrasound (QUS) measurements. Compared to DXA, QUS exposes the patient to no radiation it is non-invasive simple and easier to use, significantly less expensive and has the potential to provide a better indication of strength. QUS is particularly attractive because of its potential for widespread screening for osteoporosis. A significant number of studies have been conducted on the potential of ultrasound for evaluation of bone density and quality. Trebacz and Natali [6] report on the many studies indicating strong correlation between ultrasound and the amount of cancellous bone. Studies also reported ultrasound as good predictor of mechanical properties (shear modulus, and bone density) of cancellous bone samples is summarized [7,8]. In a recent paper, computer simulations of an adaptive damage-bone remodeling were carried out [9] to study the day-to-day activities-induced bone micro-damage. Unfortunately, the level of understanding of ultrasound interaction with bone is the primary obstacle to the advancement of its application to osteoporosis and of it becoming a fully developed diagnostic technique. Theory and experience from previous applications of QUS provide some indication to two main QUS variables, broadband ultrasound attenuation and velocity of sound, which are influenced by the mechanical properties of bone which in turn are determined by bone’s material namely density, elasticity [10]. Despite this knowledge, the complexity of sound–bone interaction has prevented the development of useful theoretical models such that only few that can robustly explain the ultrasound propagation in bone exist [5]. The mechanism(s) of how ultrasound physically interacts with bone is still unclear. Lin et al [5] report on the theoretical studies and their limitations, including the use of Biot theory and a stratified model to simulate the ultrasound propagation in trabecular bone. Consequently some theoretical and many in vitro studies have provided descriptive information about the nature of ultrasound–bone interaction but it is still unclear what role density and elasticity play in determining velocity. In this study a simple mathematical model (Sound OS) that describes the propagation of sound in bone is presented. A set of partial differential equations (PDEs) and the numerical scheme used to solve these PDEs are elaborated. Sound OS is based primarily on volumetric density and elasticity and focuses on
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
441
the propagation of velocity in the porous bone material. Therefore this study aims at investigating: (1) the sensitivity of velocity of sound on the density and elasticity of trabecular bone, and (2) the influence of bone porosity on velocity of sound.
2. Mathematical model 2.1. Problem description Trabecular architecture is a three-dimensional complex solid matrix composed of interconnected trabeculae (struts) and marrow filling the spaces between them. Based on this setup, consider an element of a trabecular bone between two pressure transducers where one is an emitter while the other is a receiver. The pressure transducer actually is a mechanical disturbance from equilibrium of the medium in front of it. In a free traveling wave one part of the medium disturbs an adjacent part, thereby imparting energy to it. This portion of the medium, in turn, disturbs another part, thereby causing flow energy in a given direction away from the source or in this case the emitting transducer. In other words, wave propagation is the result of kinetic energy at one point being transferred into potential energy at an adjacent point and vice versa. When the ultrasound enters this complex architecture, acoustic energy is partially reflected and partially transmitted through the marrow–trabeculae interface. The rate of travel of this disturbance or velocity propagation is determined by the constants of the medium. Ultrasound is propagated through a medium by means of wave motion. The medium as a whole does not move but only the disturbance. The wave can be viewed as excess pressure moving through the bone architecture at velocity ‘c’. When an ultrasound wave traverses a section of a trabecular bone its properties are affected by two primary parameters, namely the elasticity of the medium and the connectivity of the bone. The celerity or speed of the sound wave is directly proportional to the square root of the elasticity of the medium. If the elasticity of the bone increases, then so does the speed of the sound wave and consequently the VOS. The connectivity of the bone structure influences the traversing ultrasound wave by the degree and extent of angularity of the bone structure. This bone structure contributes to the attenuation of the ultrasound wave through refraction. With increased connectivity the higher the refraction and the larger the attenuation of the ultrasound wave. There are other parameters that influence the attenuation of the ultrasound wave such as absorption and the transfer of the ultrasound wave into heat added to the medium, however, they will not be considered in this study since their influence is assumed to be relatively negligible. 2.2. System of equations The propagation of ultrasound in trabecular bone can be formulated using the standard one-dimensional partial differential wave equation. This equation has been used successfully to simulate the propagation of sound in various mediums. In the present case, the partial differential wave equation is modified and reformulated to account for the parameters pertinent to our problem. These parameters are: (a) Porosity of the bone which could be indicative of the extent of bone trabeculae connectivity, (b) modulus of elasticity of the bone and (c) volumetric density of the bone. The propagation of sound in any medium can be described by a partial differential equation which governs the process. The second-order hyperbolic one-dimensional partial differential wave equation,
442
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
describing the propagation of sound waves traveling at a wave speed c is given by j2 u jt 2
= c2
j2 u jx 2
,
(1)
where u is the velocity of sound in m/s, t is time in seconds, and x is distance in meters. This equation is known as the classical wave equation. A first-order wave equation which has the properties of Eq. (1) is given by ju jt
= −c
ju jx
for c > 0.
(2)
Considering trabecular bone structure to be composed of a solid and fluid matrix, the first- and secondorder wave equations (Eqs. (2) and (1)) are modified to include the porosity of the matrix as follows: ju jt
= cf
j2 u jt 2
ju jx
= cf2
+ cs
j2 u jx 2
ju jx
+ cs2
, j2 u jx 2
(3) .
(4)
The standard celerity of the sound propagating in the fluid and bone components of the bone matrix can be expressed by the following standard established expressions, respectively: Ef , (5) cf = f
cs =
Es , s (1 − )
(6)
where Ef and Es are the modulus of elasticity for the fluid and solid parts, respectively, f and s are the density of the fluid and solid parts, and is the porosity of the bone. 2.3. Numerical scheme To solve Eqs. (3) and (4), the Lax–Wendroff finite difference scheme is used. This scheme derives the finite difference formulation from the Taylor-series expansion as follows: un+1 = uni + ut t + 21 utt (t 2 ) + 0[(t)3 ]. i
(7)
Using the wave equations (Eqs. (3) and (4)): utt = uxx (cf2 + cs2 ),
(8)
ut = −ux (cf + cs )
(9)
and substituting for utt and ut into Eq. (7) yields un+1 = uni − (cf + cs )(t)ux + 21 (cf2 + cs2 )(t 2 )uxx . i
(10)
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
443
Finally replacing ux and uxx with second-order accurate finite difference expressions gives the discretized form of equation governing the propagation of sound un+1 = uni − i
(c2 + cs2 ) (t 2 ) n (cf + cs ) (t) n (u − 2uni + uni−1 ), (ui+1 − uni−1 ) + f 2 x 2 x 2 i+1
(11)
where x is the discretized grid length and t is the simulation time step. This explicit one-step scheme is second-order accurate with a truncation error of 0[(x)2 , (t)2 ] and is stable whenever is less than or equal to 1, where =
c t x
(12)
is the courant number and c is the celerity (the speed of the sound wave front) of the propagating wave averaged for both bone and fluid matrices. This means that for any value of c across any x in the bone–fluid matrix, should not exceed 1.
3. Analysis and discussion of simulation results Eqs. (11) and (12) were programmed to solve for the time varying velocity wave front as it traverses through the bone marrow–trabeculae matrix as shown in Fig. 1. Computer simulations were done in two stages. The first stage included many simulations to analyze the sensitivity of the VOS to density, elasticity and porosity of trabecular bone. The second stage entailed comparing the present model calculations with measurements of cortical bone from sheep femur and computations from another model. To perform the simulations in stage 1, the width of the trabecular bone was specified at 4 cm to represent that of the calcaneus. The elasticity for both trabecular bone and the marrow was imposed at 2 100 000 dyn/m2 . The density of bone to fluid ratio was set at 2 to 1. The 4 cm bone was discretized into 100 smaller strips and the simulation time step was calculated by Eq. (11) to ensure computational stability. This setup gives a x equal to 0.0004 m and a t equal to 5.17 × 10−8 s. Many simulations were performed to study the interaction of density, elasticity and porosity with the average velocity of sound being the outcome variable. Using the mathematical model, the average sound
Width
Marrow-Trabeculae Matrix, ρ, ε, & porosity
Incoming Velocity of Sound, VOS
∆x
Average Velocity
Fig. 1. Numerical modeling setup.
444
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447 4000
Velocity (m/s)
3500
E2 E1.6
3000
E1.2 2500 E0.8 2000 E0.5 1500 1.4
1.5
1.6
1.7
1.8 1.9 2 Density (g/cc)
2.1
2.2
2.3
Fig. 2. Velocity–density–elasticity curves. 4000
Velocity (m/s)
3500 3000 2500 2000 1500 0.5
0.7
0.9
1.1
1.3 1.5 Elasticity
1.7
1.9
Fig. 3. Velocity–elasticity–density curves.
velocity across the 4 cm trabecular bone was calculated with changing bone density, elasticity and porosity as shown in Fig. 2 (variation of the VOS as a function of density for different elasticity values), 3 (variation of VOS as a function of elasticity for different density values) and 4 (relationship between VOS, porosity and density). Since marrow properties are close to that of water, a density of 1000 kg/m3 was used. A density of 2000 kg/m3 and an incoming wave velocity of 3680 m/s were used. In all figures, the time step was the same. This t value was found to ensure stability of calculations across all simulations. Throughout the simulations, all parameters and conditions were kept the same with the ratios of density, elasticity and porosity changing. The average velocity VOS was calculated by averaging the celerity of the calculated sound wave front at every x across the 4 cm bone thickness. In Figs. 2, 3 and 4, E represent elasticity, D represents density, all density scales are that of bone density and velocity represents the average velocity of the sound wave in the marrow–trabeculae matrix. The VOS as a function of bone density and elasticity is given in Figs. 2 and 3. Simulation results show that the relationship between VOS and density for the setup described above is linear such that a 47% change in bone density results in a 6% change in VOS (see Fig. 2). On the other hand, Fig. 3 shows that
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
445
Velocity (m/s)
5500 5000 D1.6 D1.8
4500
D2.0 D2.2
4000 3500 3000 0
0.2
0.4 0.6 Porosity
0.8
1
Fig. 4. Velocity–porosity–density curves.
the relationship between VOS and elasticity is non-linear and a 47% change in elasticity results in a 17% change in VOS. It is evident therefore that the sensitivity of VOS on bone elasticity is three times more than that of density. Fig. 4 shows the model simulation VOS results when porosity of the bone matrix is varied from 0.1 to 0.9 for different bone densities but constant elasticity. This figure depicts a non-linear relationship between VOS and porosity and may imply the as porosity increase, bone density has less effect on VOS. This effect is reduced by less than half for a 27% change in density from porosity of 0.1 to 0.9. It is clear from the simulation results that bone elasticity and porosity play a more significant role in bone strength and that role is complex by nature due to the non-linearity of the problem. The second stage of the analysis included the use of the present model to calculate the VOS obtained in the lab for sheep femoral condyle [10,5]. The available lab data included a set of VOS-porosity and another
Table 1 Comparison of model simulation results with sheep femoral condyle VOS data VOS (m/s)
Density (kg/m3 )
Porosity (%) PM
SH
PM
Elasticity (dyn/m2 ) SH
PM
SH
Measured sheep femoral condyle porosity and estimated density 1750 60 60 750e 1900 53 53 900e 2100 38 38 1100e
NR NR NR
3.67 5.00 6.50
NR NR NR
Measured sheep femoral condyle density and estimated porosity 1720 65e NR 710 1900 48e NR 920 2100 38e NR 1100
710 920 1100
3.45 5.05 6.50
NR NR NR
PM, present model calculations; SH, measured values of trabecular bone from sheep femoral condyle; NR, not reported; e, estimated.
446
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
of VOS-density values, where density–elasticity and porosity–elasticity were not available, respectively. Consequently, density was estimated for the first set using the second set and porosity was estimated for the second set using the first set (Table 1). Using the present model we were able to set the values porosity and density and perform multiple simulations by varying the elasticity of the bone matrix until the VOS measured in the lab is obtained accurately by the model and within 1% difference. It is evident from these results that the model can be effectively used to evaluate the elasticity of trabecular bone. It can also be used to evaluate the combination of any of the three variables to obtain an optimum and practical solution that matches the measured VOS. With today’s ultrasound technology for bone densitometry, this approach can be very useful to evaluate the three most important physical properties of bone for osteoporosis diagnosis and prevention. The importance of these three bone properties is due to the fact that they relate to bone strength which is directly associated with bone fracture, the primary outcome of osteoporosis.
4. Conclusion The development of the proposed mathematical model could become a useful tool in studying the interaction between ultrasound and bone architecture. The model could be used to investigate the influence of various parameters such as elasticity, density and porosity of bone on the propagation sound. The present mathematical model can be incorporated into clinical practice not only for diagnosis of osteoporosis and prevention of fractures but also for followup during treatment. When a patient is under treatment for osteoporosis, the bone has been found to increase in strength, however, this is not because of the re-establishment of the connective bone structure but due to the buildup of hydroxyapatite (calcium matrix) around the pores. This buildup can be viewed as reinforcement and not as optimal modeling of bone structure. Although treatment for osteoporosis has shown to reduce fracture risk, this might not necessarily be due only to treatment. Other factors such as lifestyle changes that may insue after a diagnosis of osteoporosis may also play a role in risk reduction by avoiding exposure to risk. In effect, the elasticity of healthy bone is not enhanced during treatment. On the one hand, the deposited bone in the bone due to the treatment has very little elasticity properties and possibly makes bone less ductile. Normal bone on the other hand, has much better elastic properties and behaves comparatively less brittle and more ductile. Much work still remains to elucidate those fundamental factors that contribute to fracture reduction after treatment.
References [1] C. Wu, H. Njeh, S. Zhao, P. Augat, D. Newitt, T. Link, Y. Lu, S. Majumdar, K.H. Genant, Ultrasound velocity of trabecular cubes reflects mainly bone density and elasticity, Calcif. Tissue Int. 64 (18) (1999) 18–23. [2] B.D. Burr, Bone material properties and mineral matrix contributions to fracture risk or age in women and men, J. Musculoskeletal Neuron Interact. 2 (3) (2002) 201–204. [3] S. Prins, H. Jorgensen, L. Jorgensen, The role of quantitative ultrasound in the assessment of bone: a review, Clin. Physiol. 18 (1) (1998) 3–17. [4] L. Mosekilde, S. Bentzen, G. Ortoft, J. Jorgensen, The predictive value of quantitative computed tomography for vertebral body compressive strength and ash density, Bone 10 (1989) 465–470.
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
447
[5] W. Lin, Q. Yi-Xian, C. Rubin, Ultrasound wave propagation in the trabecular bone predicted by the stratified model. Ann. Biomed. Eng. 29 (2001) 781–790. [6] H. Trebacz, A. Natali, Ultrasound velocity and attenuation in cancellous bone samples from lumbar vertebra and calcaneus. Osteoporosis Int. 9 (1999) 99–105. [7] R.B. Ashman, J.D. Corin, C.H. Turner, Elastic properties of cancellous bone: measurement by an ultrasonic technique, J. Biomech. 20 (1987) 979–986. [8] C.M. Langton, C.F. Njeh, R. Hodgkinson, J.D. Currey, Prediction of mechanical properties of the human calcaneus by broadband ultrasonic attenuation, Bone 18 (1996) 495–503. [9] S. Ramtani, M.J. Garcia, M. Doblare, Computer simulation of an adaptive damage-bone remodeling law applied to three unit-bone bars structure, Comput. Biol. Med. 34 (2004) 259–273. [10] S. Han, J. Rho, J. Medige, I. Ziv, Ultrasound velocity and broadband attenuation over a wide range of bone mineral density, Osteoporosis Int. 6 (1996) 291–296. Dr. Raafat George Saadé is an assistant professor at the Department of Decision Sciences and Management Information Systems, John Molson School of Business, Concordia University, Montréal, Québec, Canada. He completed his Ph.D. in 1995 and was subsequently awarded the National Research Council postdoctoral fellowship for 2 years which he completed at McGill University. Dr. Saadé has over 16 years of industrial experience and consulting in the fields of engineering, training and healthcare. Most of his academic and industrial works entail the development and use of various information systems. During the past 9 years, his work focused on the development and use of intelligent elearning/etraining and ehealth systems. Dr. Saadé’s work has resulted in numerous publications in peer refereed journals and conferences including the Journal of Information Technology in Education, Journal of Information Systems in Education, Information and Management, Expert Systems with Applications, and Computers in Human Behavior. Dr. Saadé is a founder of non-profit organization (called Viéquilibré/BalancedLife) targeting senior citizens, ‘Balanced life’ aims at screening for malnutrition. He is currently involved in a number of projects, some of the more interesting ones entail: the application of vigotsky’s triangle for inquiry-based elearning, the use of item-response theory for learning and assessment and the development of a Markovitch Genetic Algorithm hybrid system for fund management. Dr. George Tsoukas Medical Degree—McGill University (1968). In 1975 he was granted a specialization in Internal Medicine and Endocrinology—McGill. For over 10 years, he has been an associate physician at the McGill University Medical Center and is currently conducting clinical research on bone diseases. Dr. Tsoukas was also chief examiner for the Quebec College of Endocrinologists (1980–1986). Dr. Tsoukas is an educator who, has produced medical CD-ROMs, hosted a TV program explaining medical conditions and has lectured to other doctors on behalf of major pharmaceutical companies. Dr. John Caminis is a director in Clinical Development and Medical Affairs in the Bone Metabolism Section of Novartis Pharmaceuticals Corporation, East Hanover, New Jersey USA. He completed his M.D. degree in 1990 at the University of Athens in Greece. He entered a Clinical Research Fellowship at the McGill Bone Centre, Montreal General Hospital with a main emphasis on Genetics of Metabolic Bone Diseases. He stayed on as an associate researcher in the same group till 1998. Dr. Caminis then was recruited by Roche laboratories, Nutley, New Jersey, USA as a Medical Science Leader for a number of products especially those in the metabolic bone area. In May 2001 Dr. Caminis was invited to join his current company with a primary focus on Phase III drug development in metabolic bone diseases. He is also a collaborator of one of the bone team projects associated with the National Space Biomedical Researcher Institute, whose mandate is to develop countermeasures for the prevention of bone loss in microgravity. Dr. Caminis has co-authored a number of conference papers and publications he is a prolific speaker and a researcher of note.
Understanding velocity of sound in trabecular bone via computer simulations Raafat George Saadéa,∗ , George Tsoukasb , John Caminisc a Department of Decision Sciences and MIS, John Molson School of Business, Concordia University, 1455 De Maisonneuve
Building, Montreal, Queb., Canada H3G 1M8 b Center for Bone Metabolism, Faculty of Medicine, McGill University, Canada c Global Clinical Development & Medical Affairs, East Hanover, NJ 07936-1080, USA
Received 20 December 2004; received in revised form 28 February 2005; accepted 14 March 2005
Abstract Osteoporosis is a condition characterized by low bone mass and micro-architectural deterioration leading to non-traumatic fractures of the skeleton. It is a potentially debilitating condition especially for senior citizens. There is growing interest for quantitative ultrasound to measure bone mineral density. However, understanding of ultrasound–bone interaction is limited. Simulating ultrasound propagation of bone can help us better understand ultrasound–bone interaction and provide insight into bone architecture. In this study a mathematical model for the propagation of sound in bone is presented. Results demonstrate the suitability of the proposed modeling approach and the model’s capability to reproduce conditions in the lab. 䉷 2005 Published by Elsevier Ltd. Keywords: Artificial intelligence; Clinical; Disease management; Osteoporosis
1. Introduction Osteoporosis is a condition characterized by low bone mass and micro-architectural deterioration leading to non-traumatic fractures of the skeleton. It is a potentially debilitating condition especially for senior citizens and carries a direct cost of over 15 billion dollars a year in medical expenses, in the US alone [5]. Early diagnosis of osteoporosis is critical, if the risk of fracture is to be reduced by optimizing ∗ Corresponding author. Tel.: +1 450 682 5563; fax: +1 514 848 2424x2988.
E-mail address: [email protected] (R.G. Saadé). 0010-4825/$ - see front matter 䉷 2005 Published by Elsevier Ltd. doi:10.1016/j.compbiomed.2005.03.008
440
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
prevention and treatment strategies. Fracture is the principle outcome of osteoporosis, and it is essential to then to understand the determinants of fracture risk. One major determinant is intrinsic bone strength. Bone strength and its ability to resist fracture, is influenced by bone structure which includes bone size, shape, and micro-architectural characteristics [1]. Subsequently, and based on clinical observations; material properties of bone, namely density and elasticity, play an important role in risk fracture making them important surrogates for risk assessment. Bone from older individuals demonstrates decreased “tensile plastic deformation” (a result of compromised elastic properties of bone) that are attributed to molecular changes in either organic or the inorganic bone matrix [2]. This age-related reduction in the ability of bone to absorb and redistribute energy resulting from a fall may provide important information in assessing risk or diagnosing osteoporosis. Hence, the loss of energy absorption capability may be a significant factor increasing the risk of fracture in older women with low bone mass. Bone mineral density (BMD) is the most common and accepted osteoporosis diagnostic index. BMD is measured by dual energy X-ray absorptiometry, or DXA [3], and has been reported to explain about 60–75% of the variance in strength, and the remaining variance could be due to other factors such as accumulated fatigue damage, inferior bone micro-architecture and the state of bone remodeling [4]. However, the use of DXA has noted limitations: perception of excess patient exposed to radiation, expensive equipment, and its measurement of only apparent density, which is an image processing of different pixel densities, limiting its ability to proportionally correlate with strength [5]. Due to the inherent limitation of DXA, there is a growing interest in the use of quantitative ultrasound (QUS) measurements. Compared to DXA, QUS exposes the patient to no radiation it is non-invasive simple and easier to use, significantly less expensive and has the potential to provide a better indication of strength. QUS is particularly attractive because of its potential for widespread screening for osteoporosis. A significant number of studies have been conducted on the potential of ultrasound for evaluation of bone density and quality. Trebacz and Natali [6] report on the many studies indicating strong correlation between ultrasound and the amount of cancellous bone. Studies also reported ultrasound as good predictor of mechanical properties (shear modulus, and bone density) of cancellous bone samples is summarized [7,8]. In a recent paper, computer simulations of an adaptive damage-bone remodeling were carried out [9] to study the day-to-day activities-induced bone micro-damage. Unfortunately, the level of understanding of ultrasound interaction with bone is the primary obstacle to the advancement of its application to osteoporosis and of it becoming a fully developed diagnostic technique. Theory and experience from previous applications of QUS provide some indication to two main QUS variables, broadband ultrasound attenuation and velocity of sound, which are influenced by the mechanical properties of bone which in turn are determined by bone’s material namely density, elasticity [10]. Despite this knowledge, the complexity of sound–bone interaction has prevented the development of useful theoretical models such that only few that can robustly explain the ultrasound propagation in bone exist [5]. The mechanism(s) of how ultrasound physically interacts with bone is still unclear. Lin et al [5] report on the theoretical studies and their limitations, including the use of Biot theory and a stratified model to simulate the ultrasound propagation in trabecular bone. Consequently some theoretical and many in vitro studies have provided descriptive information about the nature of ultrasound–bone interaction but it is still unclear what role density and elasticity play in determining velocity. In this study a simple mathematical model (Sound OS) that describes the propagation of sound in bone is presented. A set of partial differential equations (PDEs) and the numerical scheme used to solve these PDEs are elaborated. Sound OS is based primarily on volumetric density and elasticity and focuses on
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
441
the propagation of velocity in the porous bone material. Therefore this study aims at investigating: (1) the sensitivity of velocity of sound on the density and elasticity of trabecular bone, and (2) the influence of bone porosity on velocity of sound.
2. Mathematical model 2.1. Problem description Trabecular architecture is a three-dimensional complex solid matrix composed of interconnected trabeculae (struts) and marrow filling the spaces between them. Based on this setup, consider an element of a trabecular bone between two pressure transducers where one is an emitter while the other is a receiver. The pressure transducer actually is a mechanical disturbance from equilibrium of the medium in front of it. In a free traveling wave one part of the medium disturbs an adjacent part, thereby imparting energy to it. This portion of the medium, in turn, disturbs another part, thereby causing flow energy in a given direction away from the source or in this case the emitting transducer. In other words, wave propagation is the result of kinetic energy at one point being transferred into potential energy at an adjacent point and vice versa. When the ultrasound enters this complex architecture, acoustic energy is partially reflected and partially transmitted through the marrow–trabeculae interface. The rate of travel of this disturbance or velocity propagation is determined by the constants of the medium. Ultrasound is propagated through a medium by means of wave motion. The medium as a whole does not move but only the disturbance. The wave can be viewed as excess pressure moving through the bone architecture at velocity ‘c’. When an ultrasound wave traverses a section of a trabecular bone its properties are affected by two primary parameters, namely the elasticity of the medium and the connectivity of the bone. The celerity or speed of the sound wave is directly proportional to the square root of the elasticity of the medium. If the elasticity of the bone increases, then so does the speed of the sound wave and consequently the VOS. The connectivity of the bone structure influences the traversing ultrasound wave by the degree and extent of angularity of the bone structure. This bone structure contributes to the attenuation of the ultrasound wave through refraction. With increased connectivity the higher the refraction and the larger the attenuation of the ultrasound wave. There are other parameters that influence the attenuation of the ultrasound wave such as absorption and the transfer of the ultrasound wave into heat added to the medium, however, they will not be considered in this study since their influence is assumed to be relatively negligible. 2.2. System of equations The propagation of ultrasound in trabecular bone can be formulated using the standard one-dimensional partial differential wave equation. This equation has been used successfully to simulate the propagation of sound in various mediums. In the present case, the partial differential wave equation is modified and reformulated to account for the parameters pertinent to our problem. These parameters are: (a) Porosity of the bone which could be indicative of the extent of bone trabeculae connectivity, (b) modulus of elasticity of the bone and (c) volumetric density of the bone. The propagation of sound in any medium can be described by a partial differential equation which governs the process. The second-order hyperbolic one-dimensional partial differential wave equation,
442
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
describing the propagation of sound waves traveling at a wave speed c is given by j2 u jt 2
= c2
j2 u jx 2
,
(1)
where u is the velocity of sound in m/s, t is time in seconds, and x is distance in meters. This equation is known as the classical wave equation. A first-order wave equation which has the properties of Eq. (1) is given by ju jt
= −c
ju jx
for c > 0.
(2)
Considering trabecular bone structure to be composed of a solid and fluid matrix, the first- and secondorder wave equations (Eqs. (2) and (1)) are modified to include the porosity of the matrix as follows: ju jt
= cf
j2 u jt 2
ju jx
= cf2
+ cs
j2 u jx 2
ju jx
+ cs2
, j2 u jx 2
(3) .
(4)
The standard celerity of the sound propagating in the fluid and bone components of the bone matrix can be expressed by the following standard established expressions, respectively: Ef , (5) cf = f
cs =
Es , s (1 − )
(6)
where Ef and Es are the modulus of elasticity for the fluid and solid parts, respectively, f and s are the density of the fluid and solid parts, and is the porosity of the bone. 2.3. Numerical scheme To solve Eqs. (3) and (4), the Lax–Wendroff finite difference scheme is used. This scheme derives the finite difference formulation from the Taylor-series expansion as follows: un+1 = uni + ut t + 21 utt (t 2 ) + 0[(t)3 ]. i
(7)
Using the wave equations (Eqs. (3) and (4)): utt = uxx (cf2 + cs2 ),
(8)
ut = −ux (cf + cs )
(9)
and substituting for utt and ut into Eq. (7) yields un+1 = uni − (cf + cs )(t)ux + 21 (cf2 + cs2 )(t 2 )uxx . i
(10)
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
443
Finally replacing ux and uxx with second-order accurate finite difference expressions gives the discretized form of equation governing the propagation of sound un+1 = uni − i
(c2 + cs2 ) (t 2 ) n (cf + cs ) (t) n (u − 2uni + uni−1 ), (ui+1 − uni−1 ) + f 2 x 2 x 2 i+1
(11)
where x is the discretized grid length and t is the simulation time step. This explicit one-step scheme is second-order accurate with a truncation error of 0[(x)2 , (t)2 ] and is stable whenever is less than or equal to 1, where =
c t x
(12)
is the courant number and c is the celerity (the speed of the sound wave front) of the propagating wave averaged for both bone and fluid matrices. This means that for any value of c across any x in the bone–fluid matrix, should not exceed 1.
3. Analysis and discussion of simulation results Eqs. (11) and (12) were programmed to solve for the time varying velocity wave front as it traverses through the bone marrow–trabeculae matrix as shown in Fig. 1. Computer simulations were done in two stages. The first stage included many simulations to analyze the sensitivity of the VOS to density, elasticity and porosity of trabecular bone. The second stage entailed comparing the present model calculations with measurements of cortical bone from sheep femur and computations from another model. To perform the simulations in stage 1, the width of the trabecular bone was specified at 4 cm to represent that of the calcaneus. The elasticity for both trabecular bone and the marrow was imposed at 2 100 000 dyn/m2 . The density of bone to fluid ratio was set at 2 to 1. The 4 cm bone was discretized into 100 smaller strips and the simulation time step was calculated by Eq. (11) to ensure computational stability. This setup gives a x equal to 0.0004 m and a t equal to 5.17 × 10−8 s. Many simulations were performed to study the interaction of density, elasticity and porosity with the average velocity of sound being the outcome variable. Using the mathematical model, the average sound
Width
Marrow-Trabeculae Matrix, ρ, ε, & porosity
Incoming Velocity of Sound, VOS
∆x
Average Velocity
Fig. 1. Numerical modeling setup.
444
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447 4000
Velocity (m/s)
3500
E2 E1.6
3000
E1.2 2500 E0.8 2000 E0.5 1500 1.4
1.5
1.6
1.7
1.8 1.9 2 Density (g/cc)
2.1
2.2
2.3
Fig. 2. Velocity–density–elasticity curves. 4000
Velocity (m/s)
3500 3000 2500 2000 1500 0.5
0.7
0.9
1.1
1.3 1.5 Elasticity
1.7
1.9
Fig. 3. Velocity–elasticity–density curves.
velocity across the 4 cm trabecular bone was calculated with changing bone density, elasticity and porosity as shown in Fig. 2 (variation of the VOS as a function of density for different elasticity values), 3 (variation of VOS as a function of elasticity for different density values) and 4 (relationship between VOS, porosity and density). Since marrow properties are close to that of water, a density of 1000 kg/m3 was used. A density of 2000 kg/m3 and an incoming wave velocity of 3680 m/s were used. In all figures, the time step was the same. This t value was found to ensure stability of calculations across all simulations. Throughout the simulations, all parameters and conditions were kept the same with the ratios of density, elasticity and porosity changing. The average velocity VOS was calculated by averaging the celerity of the calculated sound wave front at every x across the 4 cm bone thickness. In Figs. 2, 3 and 4, E represent elasticity, D represents density, all density scales are that of bone density and velocity represents the average velocity of the sound wave in the marrow–trabeculae matrix. The VOS as a function of bone density and elasticity is given in Figs. 2 and 3. Simulation results show that the relationship between VOS and density for the setup described above is linear such that a 47% change in bone density results in a 6% change in VOS (see Fig. 2). On the other hand, Fig. 3 shows that
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
445
Velocity (m/s)
5500 5000 D1.6 D1.8
4500
D2.0 D2.2
4000 3500 3000 0
0.2
0.4 0.6 Porosity
0.8
1
Fig. 4. Velocity–porosity–density curves.
the relationship between VOS and elasticity is non-linear and a 47% change in elasticity results in a 17% change in VOS. It is evident therefore that the sensitivity of VOS on bone elasticity is three times more than that of density. Fig. 4 shows the model simulation VOS results when porosity of the bone matrix is varied from 0.1 to 0.9 for different bone densities but constant elasticity. This figure depicts a non-linear relationship between VOS and porosity and may imply the as porosity increase, bone density has less effect on VOS. This effect is reduced by less than half for a 27% change in density from porosity of 0.1 to 0.9. It is clear from the simulation results that bone elasticity and porosity play a more significant role in bone strength and that role is complex by nature due to the non-linearity of the problem. The second stage of the analysis included the use of the present model to calculate the VOS obtained in the lab for sheep femoral condyle [10,5]. The available lab data included a set of VOS-porosity and another
Table 1 Comparison of model simulation results with sheep femoral condyle VOS data VOS (m/s)
Density (kg/m3 )
Porosity (%) PM
SH
PM
Elasticity (dyn/m2 ) SH
PM
SH
Measured sheep femoral condyle porosity and estimated density 1750 60 60 750e 1900 53 53 900e 2100 38 38 1100e
NR NR NR
3.67 5.00 6.50
NR NR NR
Measured sheep femoral condyle density and estimated porosity 1720 65e NR 710 1900 48e NR 920 2100 38e NR 1100
710 920 1100
3.45 5.05 6.50
NR NR NR
PM, present model calculations; SH, measured values of trabecular bone from sheep femoral condyle; NR, not reported; e, estimated.
446
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
of VOS-density values, where density–elasticity and porosity–elasticity were not available, respectively. Consequently, density was estimated for the first set using the second set and porosity was estimated for the second set using the first set (Table 1). Using the present model we were able to set the values porosity and density and perform multiple simulations by varying the elasticity of the bone matrix until the VOS measured in the lab is obtained accurately by the model and within 1% difference. It is evident from these results that the model can be effectively used to evaluate the elasticity of trabecular bone. It can also be used to evaluate the combination of any of the three variables to obtain an optimum and practical solution that matches the measured VOS. With today’s ultrasound technology for bone densitometry, this approach can be very useful to evaluate the three most important physical properties of bone for osteoporosis diagnosis and prevention. The importance of these three bone properties is due to the fact that they relate to bone strength which is directly associated with bone fracture, the primary outcome of osteoporosis.
4. Conclusion The development of the proposed mathematical model could become a useful tool in studying the interaction between ultrasound and bone architecture. The model could be used to investigate the influence of various parameters such as elasticity, density and porosity of bone on the propagation sound. The present mathematical model can be incorporated into clinical practice not only for diagnosis of osteoporosis and prevention of fractures but also for followup during treatment. When a patient is under treatment for osteoporosis, the bone has been found to increase in strength, however, this is not because of the re-establishment of the connective bone structure but due to the buildup of hydroxyapatite (calcium matrix) around the pores. This buildup can be viewed as reinforcement and not as optimal modeling of bone structure. Although treatment for osteoporosis has shown to reduce fracture risk, this might not necessarily be due only to treatment. Other factors such as lifestyle changes that may insue after a diagnosis of osteoporosis may also play a role in risk reduction by avoiding exposure to risk. In effect, the elasticity of healthy bone is not enhanced during treatment. On the one hand, the deposited bone in the bone due to the treatment has very little elasticity properties and possibly makes bone less ductile. Normal bone on the other hand, has much better elastic properties and behaves comparatively less brittle and more ductile. Much work still remains to elucidate those fundamental factors that contribute to fracture reduction after treatment.
References [1] C. Wu, H. Njeh, S. Zhao, P. Augat, D. Newitt, T. Link, Y. Lu, S. Majumdar, K.H. Genant, Ultrasound velocity of trabecular cubes reflects mainly bone density and elasticity, Calcif. Tissue Int. 64 (18) (1999) 18–23. [2] B.D. Burr, Bone material properties and mineral matrix contributions to fracture risk or age in women and men, J. Musculoskeletal Neuron Interact. 2 (3) (2002) 201–204. [3] S. Prins, H. Jorgensen, L. Jorgensen, The role of quantitative ultrasound in the assessment of bone: a review, Clin. Physiol. 18 (1) (1998) 3–17. [4] L. Mosekilde, S. Bentzen, G. Ortoft, J. Jorgensen, The predictive value of quantitative computed tomography for vertebral body compressive strength and ash density, Bone 10 (1989) 465–470.
R.G. Saadé et al. / Computers in Biology and Medicine 36 (2006) 439 – 447
447
[5] W. Lin, Q. Yi-Xian, C. Rubin, Ultrasound wave propagation in the trabecular bone predicted by the stratified model. Ann. Biomed. Eng. 29 (2001) 781–790. [6] H. Trebacz, A. Natali, Ultrasound velocity and attenuation in cancellous bone samples from lumbar vertebra and calcaneus. Osteoporosis Int. 9 (1999) 99–105. [7] R.B. Ashman, J.D. Corin, C.H. Turner, Elastic properties of cancellous bone: measurement by an ultrasonic technique, J. Biomech. 20 (1987) 979–986. [8] C.M. Langton, C.F. Njeh, R. Hodgkinson, J.D. Currey, Prediction of mechanical properties of the human calcaneus by broadband ultrasonic attenuation, Bone 18 (1996) 495–503. [9] S. Ramtani, M.J. Garcia, M. Doblare, Computer simulation of an adaptive damage-bone remodeling law applied to three unit-bone bars structure, Comput. Biol. Med. 34 (2004) 259–273. [10] S. Han, J. Rho, J. Medige, I. Ziv, Ultrasound velocity and broadband attenuation over a wide range of bone mineral density, Osteoporosis Int. 6 (1996) 291–296. Dr. Raafat George Saadé is an assistant professor at the Department of Decision Sciences and Management Information Systems, John Molson School of Business, Concordia University, Montréal, Québec, Canada. He completed his Ph.D. in 1995 and was subsequently awarded the National Research Council postdoctoral fellowship for 2 years which he completed at McGill University. Dr. Saadé has over 16 years of industrial experience and consulting in the fields of engineering, training and healthcare. Most of his academic and industrial works entail the development and use of various information systems. During the past 9 years, his work focused on the development and use of intelligent elearning/etraining and ehealth systems. Dr. Saadé’s work has resulted in numerous publications in peer refereed journals and conferences including the Journal of Information Technology in Education, Journal of Information Systems in Education, Information and Management, Expert Systems with Applications, and Computers in Human Behavior. Dr. Saadé is a founder of non-profit organization (called Viéquilibré/BalancedLife) targeting senior citizens, ‘Balanced life’ aims at screening for malnutrition. He is currently involved in a number of projects, some of the more interesting ones entail: the application of vigotsky’s triangle for inquiry-based elearning, the use of item-response theory for learning and assessment and the development of a Markovitch Genetic Algorithm hybrid system for fund management. Dr. George Tsoukas Medical Degree—McGill University (1968). In 1975 he was granted a specialization in Internal Medicine and Endocrinology—McGill. For over 10 years, he has been an associate physician at the McGill University Medical Center and is currently conducting clinical research on bone diseases. Dr. Tsoukas was also chief examiner for the Quebec College of Endocrinologists (1980–1986). Dr. Tsoukas is an educator who, has produced medical CD-ROMs, hosted a TV program explaining medical conditions and has lectured to other doctors on behalf of major pharmaceutical companies. Dr. John Caminis is a director in Clinical Development and Medical Affairs in the Bone Metabolism Section of Novartis Pharmaceuticals Corporation, East Hanover, New Jersey USA. He completed his M.D. degree in 1990 at the University of Athens in Greece. He entered a Clinical Research Fellowship at the McGill Bone Centre, Montreal General Hospital with a main emphasis on Genetics of Metabolic Bone Diseases. He stayed on as an associate researcher in the same group till 1998. Dr. Caminis then was recruited by Roche laboratories, Nutley, New Jersey, USA as a Medical Science Leader for a number of products especially those in the metabolic bone area. In May 2001 Dr. Caminis was invited to join his current company with a primary focus on Phase III drug development in metabolic bone diseases. He is also a collaborator of one of the bone team projects associated with the National Space Biomedical Researcher Institute, whose mandate is to develop countermeasures for the prevention of bone loss in microgravity. Dr. Caminis has co-authored a number of conference papers and publications he is a prolific speaker and a researcher of note.