Ultrasound Radiation Force for the Assessment of Bone Fracture Healing in Children: An In Vivo Pilot Study

Vibrational characteristics of bone are directly dependent on its physical properties. In this study, a vibrational method for bone evaluation is introduced. We propose a new type of quantitative vibro-acoustic method based on the acoustic radiation force of ultrasound for bone characterization in persons with fracture. Using this method, we excited the clavicle or ulna by an ultrasound radiation force pulse which induces vibrations in the bone, resulting in an acoustic wave that is measured by a hydrophone placed on the skin. The acoustic signals were used for wave velocity estimation based on a cross-correlation technique. To further separate different vibration characteristics, we adopted a variational mode decomposition technique to decompose the received signal into an ensemble of band-limited intrinsic mode functions, allowing analysis of the acoustic signals by their constitutive components. This prospective study included 15 patients: 12 with clavicle fractures and three with ulna fractures. Contralateral intact bones were used as controls. Statistical analysis demonstrated that fractured bones can be differentiated from intact ones with a detection probability of 80%. Additionally, we introduce a “healing factor” to quantify the bone healing progress which successfully tracked the progress of healing in 80% of the clavicle fractures in the study.


Introduction
Speed of fracture healing may vary from one individual to another and among different fracture locations based on anatomic location, fracture severity, patient physiological state, and treatment plan. Immobilization is the choice treatment for a majority of fractures. Determining the optimum duration of immobilization is important, as prolonged cast immobilization has multiple drawbacks, including muscle atrophy, joint stiffness, and skin breakdown. Furthermore, each cast change exposes the patient to the risk of cast saw burns or lacerations [1]. Paradoxically, excessive immobilization can delay the healing process, since bone formation needed for healing the fracture site is stimulated by mechanical stresses through mechanosensors within the osteocytes [2]. Significant debate exists amongst musculoskeletal care providers as to identifying when a fracture is healed sufficiently to discontinue immobilization and resume activities. Currently, the assessment of bone fracture healing is performed by serial radiographic and physical examinations [3]. Studies have shown the poor sensitivity of radiographic measurement in the prediction of healing in a rabbit model of bone fracture [4][5][6].
Conventional ultrasound in humans has been reported for detecting clavicle and forearm fractures with variable sensitivity [7]. A 1998 study [8] used conventional ultrasound to assess fracture healing in a preclinical canine model and found that a hyperechoic ultrasound signal correlated 100% of the time with the presence of hard fracture callus biopsy tissue. The development of quantitative tools for bone healing monitoring may assist physicians to identify the healing endpoint, thus preventing prolongation of the treatment course or incorrect timing of the removal of an external or internal fixation device. Several noninvasive quantitative techniques have been investigated to monitor bone fracture healing [9], including mechanical vibrational analysis [10][11][12] and acoustic emission [13,14]. Quantitative ultrasound (QUS) based on axial transmission has been proposed to assess normal and fractured bone [15,16].
The QUS technique involves placing a set of linearly arranged transducers (transmitters and receivers) along the bone axis. The receivers detect the ultrasound waves that have propagated along the cortical shell, and some of the generated parameters can be estimated for tissue characterization. Two-dimensional simulations on bone-mimicking plates, as well as in vivo studies, have investigated the potential of ultrasound axial transmission to evaluate fracture healing [15][16][17]. A model for fracture healing monitoring using QUS has been validated by experimental studies on ex vivo bovine femur and simulation studies [18]. Additionally, the effect of callus mineral composition has been reported by using simulation and experimental phantom studies [19]. Overall, both conventional ultrasound and QUS have variable results due to the presence of soft tissue around the bone. The impacts of overlying soft tissue on QUS and vibrational analysis in evaluating bone quality have been reported [20,21].
The utility of ultrasound radiation force (URF) for the characterization of soft tissues has been of interest for the past two decades. Vibro-acoustography, in particular, uses URF for imaging and tissue characterization [22,23]. The focused URF of vibro-acoustography induces vibrations in the object of interest. The frequency of those vibrations can be much smaller than that of the primary ultrasound used for URF. Object vibrations produce an acoustic field that is detected by a hydrophone sensitive to the frequency band of the induced vibrations. To construct an image of the object, the URF is scanned across it, usually in a raster fashion. The resulting image provides information on the object's vibrational and mechanical characteristics. The temporal waveform of the URF depends on the particular application. Examples of URF waveforms include single frequency [22,24], multi-frequency [25], tone burst [26,27], and chirp signal [28]. Depending on the application, different types of transducers have been used to produce the URF, such as dual-element confocal [29], linear array [30], multi-row array [26], and sector array [31]. Vibro-acoustography has been used for imaging tissues and implants [32][33][34][35][36]. An overview of the application of this technique on various soft tissues has been previously reported [27]. Additionally, vibro-acoustography may be used for bone evaluation, as described in a recent study by Denis and colleagues [37]. The initial application of vibro-acoustography in the assessment of bone fracture and bone healing was presented in 2006 [38]. In that study, continuous URF was used to evaluate the resonance frequencies of ex-vivo animal bone. The premise behind their method was that bone fracture and the level of repair could be determined based on changes in the bone resonance frequency(ies).
In this paper, we propose a quantitative URF method for detecting bone fracture and monitoring bone healing in human subjects. We applied focused URF excitations at four equi-distant points on the bone of interest. Resulting acoustic signals were used for wave arrival time and velocity estimation. The arrival time of the received signal at the hydrophone from each excitation included the effect of overlaying soft tissue. However, by measuring the time difference between the arrival times of the signals originating from two excitation points at known positions, it is possible to calculate the speed in bone in the vicinity of the excitation points. The signal measured by the hydrophone can be considered as a superposition of different modes or waves. Variational mode decomposition (VMD) [39] is an efficient method for decomposing a signal into an ensemble of band-limited intrinsic mode functions, and we used it to decompose the received signal into its constitutive modes and estimate the arrival time and velocity from the decomposed modes. Our proposed method successfully distinguished between fractured and normal bone with an 80% probability and was able to monitor the healing progress in 80% of the clavicular fractures.

Human Subjects
This prospective study has been conducted under an approved Mayo Clinic IRB (institutional review board) Protocol and is HIPPA (Health Insurance Portability and Accountability Act) Compliant. The study included 15 children with an obtained IRB approved informed consent/assent form, recruited between August 2016 and October 2017. Of the total, 12 had a single clavicle fracture and three had a single ulna fracture. Median age was 12 ± 4.1 years. There were 10 males (median age 11 ± 3.9 years) and five females (median age 8 ± 4.6). Signed written consent was obtained from the parents/guardians of each subject. Participants who could attend multiple visits post-fracture were considered for assessing bone healing and included 13 subjects with a clavicle fracture and three with an ulna fracture. The median times between the fracture date and the dates of the first, second, and third visits were 8 ± 2.5 days, 30 ± 8.7 days, and 45 ± 16.3 days, respectively. Clinical reports, including X-ray results and other clinical indicators reported by the physician, were used for the definitive assessment of fracture healing.

Remote Acoustic Measurements Setup
A fully programmable ultrasound platform (Verasonics V-1, Kirkland, WA, USA) equipped with a linear array transducer (L7-4) was used to guide and generate the URF excitation. The URF signal consists of a five-cycle tone burst at a 5 MHz center frequency. The magnitude radiation force exerted by an acoustic plane wave in a 5 MHz center frequency transducer depends on the absorption coefficient and average acoustic intensity wave. The absorption coefficient directly depends on the frequency and the average acoustic intensity wave directly depends on the voltage of the transducer. By increasing the frequency, the absorption coefficient is increased, while the maximum tolerable voltage of each element in the transducer is decreased. We selected the frequency and voltage in order to achieve the maximum radiation force. Among the linear array transducers (Verasonics research system), L7-4 was the best candidate for our experimental setup. A Reson TC4034-1 hydrophone (Denmark) was used to detect the radiated acoustic pressure response. The safety measurements (i.e., temperature index and mechanical index) for the setup were performed; both were within the Food and Drug Administration limit. The hydrophone output was then filtered with an analogue eight-pole band-pass Bessel filter with a center frequency of 500 kHz, digitized at a sampling frequency of 100 Mega samples per second, and recorded on a computer. The transmitting linear array and the receiving hydrophone (RX) were placed in a fixture at a distance of 0.5 cm apart. Acoustic gel was used for coupling the hydrophone/probe to the skin. The experimental setup can be seen in Figure 1.

Problem Statement
We defined two hypothesis tests: one for distinguishing between intact and fractured bones and another for identifying the progression of the healing process in fractured bone. The first hypothesis test was defined based on the velocity of the wave received by the hydrophone. The second hypothesis test was defined based on the comparison of the wave velocity distribution in the fractured bone against that of the corresponding contralateral intact bone at different visits. In the first test, the velocity of the wave received by the hydrophone was considered as a decision variable for separating fractured from intact bones, which was defined as: where H F and H I denote a set of fractured and intact bones, respectively; v is the wave velocity; and v THR is the threshold of decision based on the wave velocity [9]. This hypothesis test was defined based on the fact that we would expect the wave velocity in intact bone to be higher than the wave velocity in fractured bone. To test the hypothesis in (1), we needed to estimate the wave velocity of the received signal. Hence, in the next Sections 2.4.1 and 2.4.2, two methods are described for wave velocity estimation. For each subject, we measured the velocity of the fractured and the corresponding contralateral intact bone multiple times. Thus, for fracture detection, it was possible to compare the distributions of the velocity (or the corresponding time delay) in each pair of intact and fractured bones.

Problem Statement
We defined two hypothesis tests: one for distinguishing between intact and fractured bones and another for identifying the progression of the healing process in fractured bone. The first hypothesis test was defined based on the velocity of the wave received by the hydrophone. The second hypothesis test was defined based on the comparison of the wave velocity distribution in the fractured bone against that of the corresponding contralateral intact bone at different visits.
In the first test, the velocity of the wave received by the hydrophone was considered as a decision variable for separating fractured from intact bones, which was defined as: , where F H and I H denote a set of fractured and intact bones, respectively; v is the wave velocity; and THR v is the threshold of decision based on the wave velocity [9]. This hypothesis test was defined based on the fact that we would expect the wave velocity in intact bone to be higher than the wave velocity in fractured bone. To test the hypothesis in (1), we needed to estimate the wave velocity of the received signal. Hence, in the next Sections 2.4.1 and 2.4.2, two methods are described for wave velocity estimation. For each subject, we measured the velocity of the fractured and the corresponding contralateral intact bone multiple times. Thus, for fracture detection, it was possible The second hypothesis test was also based on comparing wave velocity distributions between fractured bone and the corresponding contralateral intact bone in each patient. To detect the progress of the healing process of each fractured bone over time, we compared the probability distribution function (PDF) of the estimated velocity in fractured bone to that of the contralateral intact bone at each visit using a rank sum test. The derived p value from the comparison of the pair of PDFs of velocity (or time delay) in patient l at visit m is denoted by p l (m). Thus, we defined a decision variable for detecting the progress of the healing process in each patient as follows: where N l is the number of visits for patient l. If p l (m) is a monotonically increasing function in terms of visit m, we detect the healing process; otherwise, we cannot consider that sufficient healing progress has occurred in patient l. Moreover, if the p value in the first visit is less than 0.05, i.e., p l (1) < 0.05, and in the last visit, p l (N l ), is sufficiently large, the decision variable of HD l = 1 is considered as healing detection; otherwise, HD l = 0. The percentage of successful healing detection was defined by the quality healing factor (QHF) as: where L is the number of patients. The advantage of comparing the PDF of velocities (or time delay) on different measurements of intact bone and bone with fracture as opposed to having only one statistic, such as the mean or median of velocity (or time delay), is that complete information on the PDF of velocity (or time delay) is maintained for comparing intact and fractured bone. Based on these two hypothesis tests and the use of the rank sum test, we examined the performance of our method for detecting fractures and fracture healing in human subjects.

Acoustic Parameters Estimation
Parameters were estimated from the temporal characteristics of the acoustic response in order to capture changes in the bone. The wave velocity and arrival time differences of the acoustic waves are of particular interest for bone mechanical assessment. For simplicity of notations, we omitted the volunteer index over all the defined parameters; however, in the results section, the distribution of each defined parameter was compared for fractured and intact bones over groups of volunteers. The detailed methods for parameter estimation used for evaluating the performance of the hypothesis tests are presented in the subsections that follow.

Wave Velocity Estimation
The velocity of the acoustic wave was estimated from the arrival time difference measurements along the bone at each visit. The probe/hydrophone assembly was positioned along the bone, such that an imaginary line between the probe and the hydrophone was parallel to the long axis of the bone. In the case of fractured bone, the fracture position was marked on the skin by the clinical staff. Using this mark as the guide, the probe/hydrophone assembly was placed over the fracture, such that the focus of the URF was on one side of the fracture while the hydrophone was on the other side of the fracture, as shown in Figure 2a. To make the measurements on the contralateral intact bone, the probe/hydrophone assembly was placed at an approximately similar position.  Multiple measurements were performed by focusing the URF at different positions on the bone. The distance between the URF focus (TX) and the receiving hydrophone (RX) is d i = d 0 + i∆d cm, where d 0 = 1.2 cm, i ∈ {0, 1, 2, 3}, and ∆d = 1.2 mm (Figure 2b). For the measured signal at distance d i , the kth measurement of mth patient visit is denoted by x (k,m) i (nT s ), where T s = 1/F s is the sampling interval and F s is the sampling frequency, as follows: (nT s ) is the n th sample at distance d i and k th measurement of m th patient visit, 1 ≤ k ≤ K, 1 ≤ n ≤ N, 1 ≤ m ≤ M. The variables K, N, and M are the number of measurements at each patient visit, the number of samples in each measurement, and the number of patient visits, respectively. We measured at K different sites on the bone. At all sites, the fracture was between the excitation point of the probe and the hydrophone.
where ∆τ , which was derived from the following optimization problem: where v max and v min are, respectively, the maximum and minimum values of the allowed velocity in the propagation medium (i.e., bone); andR where i, j ∈ {0, 1, 2, 3}, j < i, and N C are the number of samples used for the estimation of the cross correlation. Also, w[n] := w(nT s ) denotes the Tukey window (tapered cosine) in the time domain and was defined as: where 0 ≤ r ≤ 1 is the roll-off factor. To compare the bones with and without fracture, we estimated the sample mean, sample median, and sample standard deviation of the estimated time delay and bone wave velocity for each measurement and each volunteer. The mean, median, and standard deviation of ∆τ (k,m) ij at mth patient visit were derived as follows: In a similar way the mean, median, and standard deviation of v (k,m) ij at mth patient visit were obtained by: (9c)

Variational Mode Decomposition
The signal received by the hydrophone may be considered as a superposition of different components. The measured signal was decomposed into a discrete number of sub-signals (modes), u (k,m) p,i , which had specific sparsity properties related to the bandwidth of the mode. These modes are reproducing the measured signal. The bandwidth of u (k,m) p,i , using a single side-band frequency spectrum, was derived for u (k,m) p,i from the corresponding analytic signal with a Hilbert transform. By multiplying an exponential with the estimated center frequency, the mode's frequency spectrum was demodulated and transferred to the baseband. Bandwidth was estimated using a squared L 2 -norm of the gradient. Therefore, the resulting constrained variational problem for x (k,m) i (t) was obtained by minimizing the sum of the gradients of the Hilbert transform of the base band modes [39] as follows: where u is the set of center frequencies; P is the number of modes; δ(t) is the Dirac delta function; and * and 2 2 denote the convolution operation and L 2 norm, respectively. The augmented Lagrangian function for the optimization of (10) was obtained by: where α is the weight of the bandwidth constraint and λ

Accuracy of Time Delay Estimation
To assess the accuracy of the time delay, ∆τ, estimation in intact bones as the control group, we estimated the coefficient of variation (CV) of ∆τ for intact bone over all measurement sites and participant visits with different method of time delay estimation, denoted by CV total (∆τ), which was defined as: where (13b) Also, we compared CV for mth patient visit as follows: where To assess the accuracy of ∆τ (m) over different patient visits, we calculated the coefficient of variation of ∆τ (m) , denoted by CV mean (∆τ), which was defined as: where (17b)

Results
The vibro-acoustic measurements for the fractured and intact bones demonstrated significant differences in velocity and arrival times. The number of measurements, K, was equal to three in each bone. The roll-off factor of the Tukey window in the estimation of velocity and time delay, r, was 0.25. The processing window duration, T W , was considered to be 25 µs to analyze only the first part of the received signal by the hydrophone, which mainly includes the bone component wave because the wave propagates faster in bone than in soft tissue. In this section, v max = 4080 m/s and v min = 1000 m/s are considered for solving the optimization problem in (5); v max is obtained from the maximum velocity of the wave in the bone, and v min is obtained from a lower bound for the     The median of the estimated bone velocity in fractured versus normal clavicles, using ( ) m v , is plotted for comparison in Figure 5. Velocities significantly differed between groups (p < 0.05) using the rank sum test.  The median of the estimated bone velocity in fractured versus normal clavicles, using v (m) , is plotted for comparison in Figure 5. Velocities significantly differed between groups (p < 0.05) using the rank sum test. The median of the estimated bone velocity in fractured versus normal clavicles, using ( ) m v , is plotted for comparison in Figure 5. Velocities significantly differed between groups (p < 0.05) using the rank sum test. Performance of fracture detection based on a pair-wise comparison between intact and fractured bone for each patient is listed in Table 1. The best performance of fracture detection was 0.80 using the first mode of VMD as the input of the cross correlation technique.  Performance of fracture detection based on a pair-wise comparison between intact and fractured bone for each patient is listed in Table 1. The best performance of fracture detection was 0.80 using the first mode of VMD as the input of the cross correlation technique.  Table 2 shows the estimated QHF based on the time delay and velocity for different estimation methods: constrained velocity and first mode of VMD. We found that the healing process could be detected in 80% of fractured clavicle bones with this method.

Discussion
Detection of the fracture and healing process of bone is critical. In this paper, a new noninvasive and quantitative URF method was proposed for the detection of bone fracture and healing. We found that waves generated by ultrasound radiation force excitation vary in speed, depending on the healing of the broken bone to which they are applied. In particular, the velocity and arrival time of the generated waves differed significantly between intact and fractured bone. Using our method, the statistical analysis revealed an 80% detection rate in bones with fracture. In monitoring the healing process, our method, based on an estimation of the time delay, demonstrated a detection rate of 80% in the clavicle bones and around 70% of both clavicle and ulna bones. Conventional ultrasound has also been reported for the detection of clavicle and forearm fractures, with variable sensitivity [7]. The development of quantitative tools for bone healing monitoring may assist physicians in identifying the healing endpoint, thus preventing prolongation of the treatment course or incorrect timing of fixation device removal. Several noninvasive, quantitative techniques to monitor bone fracture healing have been investigated [9][10][11][12][13][14][15][16]18,19]. These studies have shown a correlation between the ultrasound wave velocity and health status of bone, therefore proposing that the ultrasound wave is an accurate method for assessing fracture.
To assess reliability in the estimation of time delay and velocity in normal bone, the coefficient of variation of time delay and velocity was calculated in the normal bone population. All measurements showed that the value was almost 20%, indicating a reliable estimation of time delay and velocity.
There are some limitations to this work. First, the bone was supposed to be an isotropic material. Consequently, the elastic coefficient was assumed to be uniform in all directions. Second, the receiving hydrophone is omnidirectional, as opposed to a directionally focused transducer. Consequently, the received pressure has contributions in all directions; however, the hydrophone is held perpendicular to the bone surface so that the displacement toward the bone (deflection) will be the main contributor to the received pressure. Third, changes in the wave velocity with respect to bone density were not investigated. Velocity is not only dependent on the elastic modulus, but is also linked to density. Fourth, using VMD, we decomposed the signal received by the hydrophone to multiple modes and analyzed the mode with the maximum power for velocity estimation. The decomposed modes, though, could have been different from the waves propagating in the bone tissue material. The results of simulation and experimental ex vivo studies on cortical bone for analyzing the lamb modes using guided waves have also been reported [41][42][43][44].
Further investigation is required to improve the classification of fractured and normal bone and to be able to adequately monitor the healing process. Additional acoustic parameters, such as attenuation, may further enhance the accuracy of classification. A larger dataset can aid in a more accurate performance in identifying fractures and detecting bone healing.

Conclusions
Waves generated by ultrasound radiation force excitation of intact bone and bone with fracture demonstrated a significant difference in the velocity and arrival time, and in the clavicle in particular. The healing process of bone with fracture was also detected by identifying the wave velocity in fractured bone over time and comparing it to corresponding intact bone. Results indicated that waves generated by ultrasound radiation force excitation vary in velocity, depending on the health of the bone to which they are applied.