Estimating the strength of bone using linear response

arXiv:cond-mat/0102217v1 [cond-mat.soft] 12 Feb 2001 Using Nonlinear Response to Estimate the Strength of an Elastic Network Gemunu H. Gunaratne Department of Physics, University of Houston, Houston, TX 77204 The Institute of Fundamental Studies, Kandy 20000, Sri Lanka Disordered networks of fragile elastic elements have been proposed as a model of inner porous regions of large bones [Gunaratne et.al., cond-mat/0009221, http://xyz.lanl.gov]. It is shown that the ratio Γ of responses of such a network to static and periodic strain can be used to estimate its ultimate (or breaking) stress. Since bone fracture in older adults results from the weakening of porous bone, we discuss the possibility of using Γ as a non-invasive diagnostic of osteoporotic bone. Osteoporosis is a major socio-economic problem in an aging population [1]. Unfortunately, therapeutic agents which can prevent and even reverse osteoporosis often induce adverse side effects [2]. Thus, non-invasive diagnostic tools to determine the necessity of therapeutic intervention are essential for effective management of osteoporosis. Bone Mineral Density (BMD), or the effective bone density is the principal such investigative tool [3]. Ultrasound transmission through bone [4] and geometrical characteristics of the inner porous region or trabecular architecture (TA) [5–7] are being studied as complementary diagnostics. tial reduction of the ultimate (or breaking) stress with decreasing BMD, and (3) a dramatic increase of bone strength following therapeutic regeneration [8]. Together they support the conjecture that elastic networks are a suitable model of mechanical properties of bone. In this Letter we use results from a numerical study of the model to introduce a possible diagnostic tool for osteoporosis. In older adults, weakening of the TA is the principal cause of increased propensity for bone fracture [4]. Analysis of models can complement mechanical studies of bone in aiding the identification of precursors of the weakening of a TA. In Ref. [8], it was proposed that a system to model mechanical properties of a TA can be obtained by adapting a disordered network of fragile elastic elements [9]. The model system includes potential energy contributions from elasticity and from changes in bond angles between adjacent springs. Furthermore, springs that are strained beyond (a predetermined value) ǫ and bond angles that change more than δ are assumed to fracture and are removed from the network. The inertia of the network is modeled by placing masses at the vertices. When in motion, each mass experiences a dissipative force proportional to its speed. Osteoporosis is modeled by random removal of a fraction ν of springs from the network, and the bone density is estimated by the fraction of remaining links. Therapeutic regeneration is introduced by strengthening springs that experience large strain (as suggested by Wolff’s law [10]). (a) (b) FIG. 1. The stress distributions on networks of size 60 × 60 with (a) ν = 10%, and (b) ν = 30% representing “healthy” and “osteoporotic” bone respectively. For clarity only the compressed bonds are shown, and darker hues represent larger stresses. Notice that the “stress backbone” of (a) is dense, while that of (b) consists of a few coherent pathways. The primary difference between “healthy” and “weak” networks is the nature of stress propagation through them. Figure 1 provides a representation of stresses in two networks subjected to uniform compression [11]. Numerical studies of the system show analogs of several mechanical properties of bone including, (1) an initially linear stress vs. strain curve that becomes nonlinear beyond the fracture of elastic elements, (2) an exponen1 required for its implementation is denoted Tν (t). ζp is chosen to be small (ζp ≪ ζ0 ), and hence Tν (t) can be assumed to be proportional to ζp . The dynamical susceptibility of the network is defined by Tν (t) = χν (t) · ζ(t). The Fourier transform of Tν (t) is given by the convolution T̂ν (ω) = χ̂ν (ω) ∗ ζ̂(ω). Since ζ̂(ω) = ζd δ(ω + Ω0 ), it follows that χ̂ν (ω) = T̂ν (ω − Ω0 )/ζd [15]. Figure 2 shows the power spectra |χ̂ν (ω)|/χ0 for values of ν of 0%, 10%, 20%, 30%, and 40%. Though both the static and dynamical susceptibilities reduce with advancing “osteoporosis” (increasing ν), χν (t) experiences a smaller reduction. This is possibly due to the formation of additional (temporary) stress pathways when the network is subjected to periodic strain. Figure 1(a) shows the effect on a “healthy” network (ν = 10%) where large stresses supporting propagation are seen to form a dense subset. In contrast, elastic elements supporting a “weak” network (ν = 30%) form a few coherent pathways (Figure 1(b)). We refer to the set active in stress propagation as the “stress backbone” of a network [12]. For a wide range of control parameters, it is seen to become finer with increasing ν. It is easy to understand how the nature of the stress backbone is related to the stability of a bone. Loss of connectivity of a healthy TA (due to trauma) will have little effect on its stability, since many alternative pathways are available for stress propagation. In contrast, fracture of a critical link (i.e., one belonging to the stress backbone) in a weak network will have a significant impact on the remaining stress pathways, likely inducing further fracture of elastic elements. Below, we discuss how these variations in stress backbones effect the nonlinear response of networks to externally applied strain. Consider an elastic network from which a fraction ν of elastic elements have been removed. It is first subjected to an adiabatically reached compression ζ0 [13]. ζ0 is chosen well below the yield point so that there is no fracture of elastic elements. Under these conditions, the stress T0 needed to maintain the compression is proportional to ζ0 [14,8], and the static susceptibility of the network is defined as χ0 = T0 /ζ0 . 12 8 4 0 0 12000 |χ (ω)| / χ ν 0 1 105 × Γ(ν) 2 FIG. 3. The relationship between the ultimate (or breaking) stress U (ν) and Γ(ν) for several elastic networks subjected to different compressions ζ0 . The symbols ‘+’ (blue) and ‘x’ (green) represent distinct networks with identical control parameters compressed by ζ0 = 4.0 and ζ0 = 8.0 respectively. The symbols ’o’ (red) and ‘⋄’ (aquamarine) represent two other network generated with different sets of control parameters, and compressed by ζ0 = 5.0 and ζ0 = 3.0 respectively. For these networks the values of U (ν) have been rescaled as described in the text. 8000 4000 0 −1 U(ν) 0 1 Variations of these power spectra can be quantified using the ratio Γ(ν) defined by Z 1 −1 Γ (ν) = |χ̂ν (ω)|dω, (1) χ0 2 (x−Ω )/Ω 0 0 FIG. 2. The power spectra of the dynamical susceptibility (with ζd = 0.001 and Ω0 = 100) of several elastic networks normalized by the static susceptibility. The curves correspond to values of ν of 0% (blue), 10% (green), 20% (red), 30% (aquamarine), and 40% (magenta). the range of the integral being a Nyquist frequency domain [16]. Figure 3 shows a remarkable relationship between the ultimate stress U (ν) and Γ(ν) for several networks. The symbols ‘+’ and ‘x’ represent two distinct networks (constructed using different random seeds) with Next, this compressed network is subjected to an additional periodic strain ζ(t) = ζp exp(iΩ0 t), and the force 2 identical control parameters [11] which have been subjected to compressions of ζ0 = 4.0 and ζ0 = 8.0 units respectively. U (Γ) is seen lie on a common bi-linear curve. The behavior of U (Γ) for a third network is included to determine the effects of increasing the range of elastic moduli of the network [17]. The ultimate stress of a network is expected to reduce when the range of spring constants is increased (say, by a factor f ). This is because of the increase of the number of weak springs in the network. To compensate for this decrease, we rescale U (ν) (heuristically) byf . Then, symbols ‘o’ representing the third network fall on the curve determined by the first pair of networks. U (Γ) will, of course, depend linearly on the mean elastic modulus, as has been confirmed. Finally, we include results from a fourth network to study variations of the fracture strain of elastic elements. A fixed-strain fracture criterion was used in the model; i.e., any spring that is strained by more than ǫ, and and bond angle that is changed by more than δ are removed from the network. This choice was motivated by observations from recent mechanical studies which show that a trabecular architecture from a given location fails at a fixed level of strain, independent of the strength or elastic modulus of the bone [18]. However, bone samples from distinct locations exhibit different levels of fragility. For example samples from proximal rat tibia, human tibia, and human lumbar spine have been shown to fracture at strain levels of approximately 5%, 1% and 7% respectively [18–20]. The first three networks represented in Figure 3 included a common fracture criterion; specifically ǫ = ǫ0 = 5% and δ = δ0 = 0.1. Since Γ(ν) is independent of ǫ and δ and U (ν) can be expected increase with them, U (Γ) will depend on ǫ and δ. The symbols ‘⋄’ in Figure 3 represent a fourth network [21], where the ultimate stress U (ν) has been rescaled by a factor (ǫ0 /ǫ1 ), ǫ1 being the new value of the fracture strain. For parameters chosen, failure of elastic elements (as opposed to bond-breaking) was the dominant (though not exclusive) mode of fracture, justifying the use of this rescaling factor [22]. Once U (ν) is rescaled, points representing all networks collapse to the same bi-linear curve [23]. The transition between the two linear segments of U (Γ) is accompanied by a qualitative change in the stress distribution on compressed networks. Points on the right segment of each curve represent “healthy” networks; i.e., those with smaller values of ν. The histograms of stresses for these networks, an example of which is shown by the solid line of Figure 4, exhibit broad peaks. The presence of such peaks is consistent with the geometry of their stress backbones, see Figure 1(a). In contrast, points Number of Springs on the left segment of U (Γ) represent “weak” networks (i.e., larger ν), whose stress histograms exhibit no broad peaks, and have exponential tails [24]. This behavior of the histograms is consistent with the presence of a sparse stress backbones. The transition between the two linear segments of U (Γ) with increasing ν coincides with the elimination of the peak in the histogram. 3 10 2 10 1 10 −0.07 −0.05 −0.03 −0.01 0.01 Stress FIG. 4. The solid line shows a histogram of stresses for a compressed network modeling a healthy TA (ν = 10%). The dashed line shows the analog for a weak TA (ν = 25%). Each histogram represents the average of 50 configurations with identical control parameters. The elimination of the broad peak of the histogram coincides with the transition between the two linear segments of U (Γ). Positive and negative values of the stress denote extended and compressed springs respectively. The ultimate stress U of a bone is the measure essential for management of osteoporosis. Unfortunately, it is not accessible in-vivo (without breaking a bone!), and surrogates such as the bone density are used to estimate U . The BMD of a patient is compared with that of a sample population to determine if and when therapeutic interventions are necessary. However U is known to depend on other factors of bone such as the architecture of its TA and the “quality” of bone material. The significant variations introduced by these factors makes it difficult to identify individuals susceptible to fracture using measurements of BMD alone [25]. These problems can be avoided if a characteristic that relates U to an intrinsic property (i.e., one that does not need to be compared to that of a population) of a bone is available. Since factors such as bone quality and architecture of the porous medium effect both U and Γ it is conceivable that the relationship between them remains unchanged between patients. Numerical analyses of our 3 model system justify this expectation when U is rescaled by a location dependent factor (quantifying the stiffness and fracture strain of trabecular elements). Vibrational analysis has been used for in-vivo studies of bone strength and protocols are available to obtain measurements needed to evaluate Γ [26,27]. Once the stiffness and fracture strain of different bone locations are tabulated, rescaling factors required for Figure 3 can be determined. Subsequently, the ratio Γ of responses of a bone to stationary and periodic strain can be used as a non-invasive diagnostic of bone strength. The author would like to thank S. R. Nagel for pointing out that nonlinear response is related to the stress backbone. He also acknowledges discussions with M. P. Marder and S. J. Wimalawansa. This research is partially funded by the National Science Foundation, the Office of Naval Research and the Texas Higher Education Coordinating Board. [12] [13] [14] [15] [16] [17] [18] [1] P. Kiberstis, O. Smith, and C. Norman, Science, 289, 1497 (2000). [2] R. S. Weinstein, J. Bone. Miner. Res., 15, 621 (2000). [3] C. E. Cann, Radiology, 140, 813 (1981). [4] C. F. Njeh, D. Hans, T. Fuerst, C.-C Glüer, and H. K. Genant, “Quantitative Ultrasound: Assessment of Osteoporosis and Bone Status,” Martin Dunitz, London (1999). [5] L. Pothuaud, C. L. Benhamou, P. Porion, E. Lespessailles, R. Harba, and P. Levitz, J. Bone. Miner. Res., 15, 691 (2000). [6] E. Legrand, D. Chappard, C. Pascaretti, M. Duquenne, S. Krebs, V. Rohmer, M-F. Basle, and M. Audran, J. Bone. Miner. Res., 15, 13 (2000). [7] S. Majumdar, D. Newitt, M. Jergas, A. Gies, E. Chiu, D. Osman, J. Keltner, J. Keyak, and H. K. Genant, Bone, 17, 417 (1995). [8] G. H. Gunaratne, K. E. Bassler, K. K. Mohanty, and S. J. Wimalawansa, “A Model for Bone Strength and Osteoporotic Fracture,” cond-mat/0009221 at http://xyz.lanl.gov. [9] J. W. Chung, A. Roos, J. Th. M. De Hosson, and E. Van der Giesses, Phys. Rev. B, 54, 15094 (1996). [10] S. C. Cowin, A. M. Sadegh, and G. M. Luo, J. Biomech. Eng., 114, 129-136 (1992). [11] These networks were obtained by randomly displacing vertices by up to 1 unit from a square lattice with sides of 10 units. The elastic and bond-bending spring constants were chosen within [0.5, 1.5] units and [2.5, 7.5] units [19] [20] [21] [22] [23] [24] [25] [26] [27] 4 respectively, while ǫ = 5% and δ = 0.1. Each mass was 1 unit, and the dissipation coefficient η = 0.1. Distinct networks with these control parameters were generated by using different random seeds. C. Moukarzel and P. M. Duxbury, Phys. Rev. Lett., 75, 4055 (1995). The stationary states of the system are evaluated using the conjugate gradient method to minimize potential energy. The dynamical properties of the network are calculated using the Bulirsch-Stoer method [16]. Y. C. Fung, “Biomechanics: Mechanical Properties of Living Tissue,” Springer-Verlag, New York (1993). Even though χν (t) and Tν (t) can depend on the driving frequency Ω0 , numerical studies show no such dependence when ω is scaled by Ω0 . W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vettering, “Numerical Recipes - The Art of Scientific Computing,” Cambridge University Press, Cambridge, 1988. For this network, elastic and bond bending spring constants were chosen within [0.3, 1.7] units and [0.4, 3.6] units respectively, while η = 0.001 and ζ0 = 5.0 units. The exponential decay of U (ν) was significantly faster for this network. H. A. Hogan, S. P. Ruhmann, and H. W. Sampson, J. Bone Miner. Res., 15, 284 (2000); C. M. Ford and T. M. Keaveny, J. Biomech., 29, 1309 (1996). L. Rohl, E. Larson, F. Linde, A. Orgaard, and J. Jorgensen, J. Biomechanics, 24, 1143 (1991). L. Mosekilde, L. Mosekilde, and C. C. Danielsen, Bone, 8, 79 (1987). The elastic and bond bending spring constants for this network varied within ranges [0.5, 1.5] and [2.5, 7.5] respectively, while η = 0.01 and ζ0 = 3.0. Fracture criteria are given by ǫ = 2% and δ = 0.04. To the best of the author’s knowledge, there has not been a study of fracture of single trabeculae or bond-angles. The dominant mode of fracture (critical strain vs. bond breaking) of TAs from bone samples is also not known. The standard error for the left and right segments were 0.17 and 0.40 respectively. In the absence of an independent estimate of standard deviations, it is not meaningful to calculate the goodness of fit for the least square estimates. See Ref. [16] for details. S. Chan and J. Macha, Phys. Rev. B, 49, 120 (1994). D. Marshall, O. Johnell, and H. Wedel, BMJ, 312, 12541259 (1996); L. J. Melton, E. A. Chrischilles, C. Cooper, A. N. Lane, and B. L. Riggs, J. Bone. Miner. Res., 7, 1005 (1992). B. Couteau, M.-C. Hobatho, R. Darmana, J.-C. Bribnola, and J.-Y. Arlaud, J. Biomech., 31, 383 (1998); G. Van der Perre, R. Van Audekerke, M. Martens, and J.-C. Mulier, J. Biomech. Eng., 105, 215 (1990). J. J. Kaufman and T. A. Einhorn, Osteoporos. Int., 8, 517 (1993).