Characterization of bone microstructure using photoacoustic spectrum analysis.

Osteoporosis is a progressive bone disease that is characterized by a decrease in bone mass and the deterioration in bone microarchitecture. This study investigates the feasibility of characterizing bone microstructure by analyzing the frequency spectrum of the photoacoustic (PA) signal from the bone. Modeling and numerical simulation of PA signal were performed on trabecular bone simulations and CT scans with different trabecular thicknesses. The resulting quasi-linear photoacoustic spectra were fittted by linear regression, from which the spectral parameter slope was quantified. The simulation based on two different models both demonstrate that bone specimens with thinner trabecular thicknesses have higher slope. Experiment on osteoporotic rat femoral heads with different mineral content was conducted. The finding from the experiment was in good agreement with the simulation, demonstrating that the frequency-domain analysis of PA signals can provide an objective assessment of bone microstructure and deterioration. Considering that PA measurement is non-ionizing, non-invasive, and has sufficient penetration in both calcified and non-calcified tissues, this new bone evaluation method based on photoacoustic spectral analysis holds potential for clinical management of osteoporosis and other bone diseases.


Introduction
Osteoporosis, widely recognized as a major health issue, is a progressive bone disease that is characterized by a decrease in bone mass and density which can lead to an increased risk of fracture. In osteoporosis, the bone mineral density (BMD) is reduced, bone microarchitecture (BMA) deteriorates, and the amount and variety of proteins in bone are altered. Todate, most of the clinically available diagnostic methods are based on the use of either X-ray or ultrasound (US) [1]. The BMD provided by dual X-ray absorptiometry (DXA) is considered as the "gold standard" for osteoporosis diagnosis. BMD, however, can only explain 60-80% of the variability in bone strength [2]. It has been demonstrated that other mechanical factors, including microarchitecture and post-yield mechanical properties that cannot be tested by DXA, are also important in determining the fracture risk of bone [3,4]. Non-ionizing and non-invasive quantitative ultrasound (QUS) technologies provide a practical and low-cost alternative for DXA. QUS assessment of bone structure and strength is mainly performed in the transmission mode and based on the measurement of two key parameters including speed of sound (SOS) and broadband ultrasonic attenuation (BUA). These parameters are strongly correlated with BMD but less reflective of the information of bone microstructure [5][6][7].
Photoacoustic (PA) imaging and sensing, an emerging hybrid technology involving both light and sound, has excellent sensitivity to the chemical and physiological information in biological samples and has been explored for potential application in the diagnosis of osteoporosis [8][9][10][11][12][13][14][15]. Most of the previous research on PA evaluation bone were focused on either the chemical components or ultrasound properties of bone. There is no thorough research yet to explore the feasibility of PA technique in characterizing trabecular bone microstructure which is highly related to bone health. Recently studies from our group and others have demonstrated that the frequency domain power distribution of radio-frequency (RF) PA signal contains information representative of histological microstructures in biological samples [16][17][18][19][20]. To quantitatively evaluate microstructures for potential tissue characterization, a new technique named photoacoustic spectral analysis (PASA), based on frequency domain analysis of PA signals, has been extended from the framework of ultrasound spectral analysis. To extract the main features of the PA signal from a target tissue in the frequency domain, the power spectrum of the radio-frequency PA signal is computed. The power spectrum, after calibration by removing the system impulse response, is fitted by linear regression which leads to three spectral parameters including midband fit, intercept, and slope. The PA spectral parameters are highly relevant to the microstructures of optical absorbers in phantoms and biological tissues [16,17]. Imaging and sensing the PA spectral parameters also show protential for assessment of cancerous tumors and liver conditions [20,21]. PASA may offer fundamental advantages for addressing a number of practical problems faced by conventional PA imaging. For example, PASA separates the effects from system components and tissue properties on image features, and delivers system-independent quantitative results. Moreover, performing linear fitting of the averaged power spectrum provides a cogent means of addressing the stochastic nature of tissue microstructures, and can lead to quantitative and repeatable measurements.
In this study, we have, for the first time, examined the feasibility of the newly developed PASA technique in characterizing the microstructures of trabecular bone, and in differentiating bone loss and preservation from normal. Both theoretical modeling and experimental measurements on well-established animal models were conducted. Compared to the other two spectral parameters (i.e. midband fit and intercept), slope is more sensitive to the tissue heterogeneity, and can better reflect the histological microstructures and the spatial distributions of the optical absorbing chemical components in porous trabecular bone. Moreover, slope is less susceptible to the variation in light fluence and, hence, is more reliable for quantitative imaging and objective tissue characterization. Therefore, in this study, the PA spectral parameter slope of each bone specimen was quantified and correlated with the BMD of the bone.

Simulation
To evaluate the feasibility of PASA in characterizing trabecular bone microstructures, we first conducted simulations on bone models. Two different models both reflecting the changes in bone microarchitectures were used. The first model, as shown in Fig. 1(A), was based on artificial bones with different levels of bone loss. The porous feature in sample 6 was generated by a computer, with the white area showing the trabeculae and the dark area showing the pores. To represent bone loss at different severities, the samples from 1 to 5 experienced different levels of erosion which gradually reduced the trabecular thickness. The mean trabecular thickness (MTT), for bone samples 1-6 were 0.24, 0.27, 0.31, 0.37, 0.44 and 0.51 mm respectively. The MTT of each sample was quantified following the indirect indices of mean trabecular plate thickness method [22]. The RF PA signal from each sample in Fig.  1(A) were simulated using the k-wave toolbox of Matlab which was developed by the Cox group [23]. To reduce the computational cost, the simulations were performed in a twodimensional space. However, the conclusions from two-dimensional and three-dimensional simulations are the same, because both follow the same relationship between the PA spectral parameters and the physical sizes of optical absorbers [17]. In simulation, the speed of sound and the density were set as 1500 m/s and 1000 kg/m 3 , and 3200 m/s, 1900 kg/m 3 , respectively, for water and trabecular bone [23]. The optical absorption coefficient and the Gruneisen parameter of trabeculae are both relevant to the absolute value of the PA signal amplitude. However, since none of them affects the shape of PA signal waveform, these two parameters do not impact the PA spectral parameter slope. Therefore, both of them were set as 1 in simulation. For each bone sample, we simulated the RF PA signals that are received at 50 different positions evenly distributed around the sample. Then the power spectra of the 50 RF PA signals were averaged. The averaging over 50 positions can further emphasize the main features of the power spectrum, and makes the linear fit and the PA spectral parameters less susceptible to the stochastic tissue microstructures. In Fig. 1(B), the normalized power spectrum of the six bone samples following different levels of bone loss are compared. As the bone loss increases (i.e. more bone erosion), the high frequency components of the PA power spectrum increase, which meets our expectation well. When trabeculae are the dominant PA sources in trabecular bone, thicker trabeculae produce PA signal with lower high frequency components. By performing a linear fit of each power spectrum, the spectral parameter slope was quantified. In Fig. 1(C), the slope as a function of the MTT, as well as the rational fitting, is presented. As expected, smaller MTT (or severer bone loss) leads to higher spectral parameter slope.
The second bone model for simulation was based on the high resolution micro-CT images acquired from three rat bones with different BMD. The details of rat models of bone loss and preservation are in section 3.2. As shown in Fig. 2(A), the micro-CT images were from three rat tibia (OVX, Sham, and OVX + ZOL). Each amplitude image was binarized, as shown in Fig. 2(B), with the white area showing the trabeculae and the dark area showing the pores. They represent different levels of MTT (0.050, 0.058, 0.061 mm, respectively). The second model for simulation offers only three different levels of trabecular thicknesses. However, compared to the first artificial bone model, the scond model, which is based on real bone specimens, is more realistic and may better reflect the spongy feature of trabecular bone. The RF PA signals from each bone sample in Fig. 2(B) were also simulated using the same Matlab k-wave toolbox [23]. With the RF PA signals received at 50 positions evenly distributed around the sample, the power spectra were averaged, as shown in Fig. 2(C), and a linear fit of the power spectra was performed to draw the quantified PA spectral parameter slope of each bone sample. Same as the finding from the first model, the simulation from this second model suggests that thinner MTT leads to higher spectral parameter slope again, as shown in Fig.  2(D).

Experiment setup
As show in the schematic of the experiment setup in Fig. 3(A), an OPO system (Vibrant B, Opotek) pumped by an Nd:YAG laser (Brilliant B, Bigsky) was used to provide laser pulses with a repetition rate of 10 Hz and a pulse width of 5.5 ns. The laser beam with 2 mm in diameter illuminating the bone surface generates PA signals which were received by a 20 MHz focused transducer (V317, Panametrics). Working at 685-nm wavelength, the light fluence on the bone surface was 19.2 mJ/cm 2 which was within the ANSI safety limit. The bone sample and the transducer were immersed in a water bath for acoustic coupling. The PA signals from the bone specimen were recorded by a digital oscilloscope (TDS 540B, Tektronics). An example PA signal from a rat femur is shown in Fig. 3(B). This A-line includes both signal from the trabecular architecture and signal from the cortical structure. In this study, we tried to avoid the strong signal from the cortical bone and focus on the trabecular bone only. To achieve this, the direction of light illumination and the direction of PA signal detection were arranged with an angle between them larger than 90 degrees. In this case, the strong PA signal from the cortical bone at the light illuminated side arrived at the transducer later than the PA signal from the trabecular part. As shown in Fig. 3(B), the trabecular signal can be differentiated from the signals generated in cortical layers based on their difference in time of flight.
The calibration measurement used for PA spectral analysis was performed on a hair fiber with a diameter of 70 μm, as shown in Fig. 3 (C). The PA signal from the hair fiber was generated using the same light illumination geometry as those for the measurement on bone specimens. The hair fiber was placed at the focal point of the transducer. With a small size comparable to the acoustic wavelength at the center frequency of the transducer, the hair fiber illuminated by a narrow laser beam formed a point source leading to a broadband PA signal, as shown in Fig. 3(D). To account for the system response, the measurement from each bone was calibrated by dividing the PA power spectrum from the bone by the PA power spectrum from the point source. The PSD of the PA signal in Fig. 3 (B) before and after calibration was shown in Fig. 3(E). After calibration for each bone sample, the linear regression was performed in the spectral range of 2.8-31.5 MHz which was determed by the −30 dB of the power spectrum from the point source, as shown in Fig. 3(D).

Animal model
In this study, well-established rat models of bone loss and preservation were employed. 3-4 month female Sprague-Dawley rats were subject to sham surgery (Sham, N = 4), ovariectomy (OVX, N = 4), or OVX plus weekly Zoledronic Acid (OVX + ZOL, N = 4, 1.6 µg/kg/wk i.p.) [24]. The OVX rat model is required by the FDA for the evaluation of agents used to treat or prevent postmenopausal osteoporosis, and the model has been validated based on the early bone turnover that produces bone loss following estrogen withdrawal. 4 weeks after surgery, rats were euthanized, and femora were dissected and subject to PA assessment. Significant reduction in bone mass in the femoral heads was evident in the OVX rats, while rats from the OVX + ZOL group provide a model of bone preservation. Using the setup shown in Fig. 3, PA measurement was performed ex vivo on the femoral head. To validate the findings from PA measurements, the bone specimens were also scanned by a micro-CT (GE Healthcare eXplore Locus RS). The microCT images, as shown in Fig. 4, verified that the MTT of OVX diminished significantly, resulting in 56.0% reduction; while for OVX + ZOL rats, the MTT increased about 18.1% over the sham controls.   , as well as the east-square linear fit, from the three bone groups (i.e. OVX, Sham, OVX + ZOL). The linear fit enabled the extraction of the spectral parameters including the slope that were further analyzed. In comparison with the sham control (normal), the power spectra from the OVX bones containing less and thinner trabeculae have stronger high frequency components; while the OVX + ZOL bones containing higher amount and larger size of trabeculae show weaker compared in Fig. 5(B), suggesting that, at 685-nm wavelength, the osteoporosis bones (OVX) have larger slope than the bones from the other two groups with higher BMD. Although with considerable overlap, the measurements from the sham and the OVX + ZOL groups also show noticeable difference. The finding from the experiment on the rat models of bone loss and preservation has a good agreement with that from the simulations.

Conclusion and discussion
The results from this study suggest the feasibility of a novel PASA technique, by quantitatively analyzing the frequency spectrum of the RF PA signal measured from the bone, in objective assessment of bone microarchitecture. The MTT, an important factor determining the bone strength, can be correlated with the spectral parameter slope quantified from the PA power spectrum. Considering the unique advantages of PA measurement, including high sensitivity, being non-ionizing and non-invasive, and sufficient penetration in both calcified and non-calcified tissues, this new PASA technique holds potential for clinical translation to clinic. However, studies involving larger number of samples from animal models and human patients are necessary to verify these findings.
In this work exploring the relationship between PA spectral parameter and bone microstructure, simulations based on k-space algorithm were conducted, which helped us to understand the findings from the experiment. k-space algorithm can deal with the complications during PA signal propagation through the bone including ultrasonic scattering, reflection, and refraction at the numerous interfaces. However, k-space algorithm does not consider shear waves, various propagating acoustic modes, frequency dependent acoustic velocity and material attenuation which all affect the propagation and detection of PA signal from the bone. Therefore, future study based on more accurate modeling is needed to fully understand how each bone parameter affects PA signal in time and frequency domain.
As another limitation of this work, neither the simulation nor the experiment has considered the soft tissue covering the bone. As a major impact, the ultrasonic attenuation in the soft tissue overlying the bone may change the PA power spectrum and the spectral parameter slope. As an example, considering that the mean thickness of the soft tissue over human calcaneus bone is 8.5 +/− 1.5 mm [25] and the ultrasonic attenuation coefficient in soft tissue is about 0.54 (dB/(MHz·cm)) [26], the PA spectral parameter slope is estimated to decrease about 0.46 +/− 0.08 (dB/MHz) as a result of the ultrasonic attenuation in overlying soft tissue. When the thickness of overlying soft tissue in known, this impact, however, can be largely compensated, especially when PASA is performed in conjunction with QUS [27].
Considering that both PASA and QUS involve broad band ultrasonic detection, PASA is a natural complement to QUS techniques which have been extensively studied and have already led to clinical instruments. In future clinical settings, PASA could be combined with QUS to realize PA-QUS dual measurement of bone. It is worth noting that the parameter slope of BUA quantified in QUS reflects the relationship between the attenuation of sound signals through bone and its frequency. The spectral parameter slope as quantified in PASA refers to the descent rate of the PA spectral magnitude with respect to the frequency. Therefore, these two shopes, although with the same unit of dB/MHz, have different physical meanings, and may have different correlations with bone microstructure. Unlike QUS, the ultrasonic signal in PASA is not produced by the transmit probe but instead by the laser light illumination in bone. The frequency components of the initial pressure waves in a trabecular bone are mainly determined by its microstructal features. Therefore, it is reasonable to expect that the PA signal from bone could carry richer and more direct information of bone microstructure. Before being received by the transducer, the PA signal from a bone would also be affected by the ultrasonic attenuation in the trabecular and cortical parts of the bone as well as in the overlying soft tissues. In bones with higher BMD, the ultrasound attenuation is stronger, especially in the high frequency range, which also affects the shape of the power spectrum of the received PA signal. Therefore, the PA spectral parameter slope is a combination of not only bone microstructures but also its ultrasonic properties. In order to separate these two aspects, the bone ultrasonic properties which could be quantified by QUS may be introduced into the procedure of PASA to achieve better accuracy for bone characterization. In comparison with DEXA, QUS has lower device and operating costs. In this proof-of-principle study, the light source for PASA was from an Nd:YAG laser pumped OPO which significantly increased the cost of this bone evaluation device. In the future, the PASA of bone can be implemented inexpensively by using high-power laser diodes which has already been adapted to biomedical PA imaging and sensing [28,29].