A non-contact galvanometer-based optical scanning system for diffuse correlation tomography was developed for monitoring bone graft healing in a murine femur model. A linear image reconstruction algorithm for diffuse correlation tomography was tested using finite-element method based simulated data and experimental data from a femur or a tube suspended in a homogeneous liquid phantom. Finally, the non-contact system was utilized to monitor in vivo blood flow changes prior to and one week after bone graft transplantation within murine femurs. Localized blood flow changes were observed in three mice, demonstrating a potential for quantification of longitudinal blood flow associated with bone graft healing.
© 2015 Optical Society of America
Annually over 500,000 and 2.2 million bone graft procedures are performed in the United States and worldwide, respectively [1, 2 ]. Although transplantation of autologous bone grafts (i.e., autografts), isolated from non-load-bearing regions of the skeleton, offer excellent healing, tissue volume is limited by donor site morbidity. Therefore, for critical-sized bone defects requiring large tissue volumes, the current clinical gold standard involves the use of processed cadaveric allografts . Allografts, however, result in 60% 10-year post-implantation failure rates due to poor integration and remodeling into host tissue, which is due in part to poor vascularization [3–5 ].
Blood flow is important for multiple phases in bone healing since blood vessels deliver oxygen, nutrients, growth factors, and circulating cells to the graft sites . However, skeletal blood flow has not been widely studied due to limitations with existing techniques. Labeled microsphere administration is generally considered the gold standard technique for blood flow quantification, but is regarded as experimentally difficult and requires sacrifice of animals . Laser Doppler imaging is a non-invasive technique to measure two-dimensional blood flow, but is limited due to shallow penetration depths . Positron emission tomography is another in vivo imaging technique that provides a three-dimensional map of absolute flow rates, but is expensive and subject to the availability of radioisotopes (e.g., 15O water) [8, 9 ]. Doppler ultrasound is limited in quantifying blood flow velocity of blood vessels smaller than 100 ~200 μm [10–12 ]. Furthermore, measurement of skeletal blood flow by Doppler ultrasound is not feasible because the bone acts as an acoustic barrier .
Diffuse correlation spectroscopy (DCS) and tomography (DCT) are deep-tissue, non-invasive techniques that quantify microvascular blood flow . They are suitable for longitudinal monitoring studies since they do not use ionizing radiation or contrast agents and are inexpensive to implement. DCS, which quantifies bulk blood flow with the simplified assumption of homogeneous media, has been applied to tumors, brain and skeletal muscle in animals and humans in vivo [15–21 ]. Furthermore, it has been validated by comparison with laser Doppler flowmetry [22, 23 ], Doppler ultrasound [17, 24–27 ], arterial-spin labeled magnetic resonance imaging [28–31 ], xenon computed tomography (CT) , and fluorescent microsphere measurements . DCT, which provides three-dimensional blood flow images, has been limited to applications in rat brain [33, 34 ], due in part to the difficulties of obtaining large spatial data sets necessary for reliable image reconstruction and good signal-to-noise ratios.
In this study, we developed a non-contact system to provide DCT of murine femurs, towards the goal of monitoring temporal and spatial blood flow changes during bone graft healing. Several types of non-contact DCS or DCT systems have been developed and utilized by others [33–36 ]. One type projected an image of a fixed distribution of source and detector fibers onto the tissue surface through a conventional camera with crossed polarizers for rat brains and murine tumors [17, 33, 34 ]. Another type featured the separation of optical paths between sources and detectors for larger source-detector separations suitable for human tissue applications  and its extension to collect large spatial data sets using a linear motorized stage . Our approach is unique in that it utilizes galvanometer-based optical scanning which enables collection of large spatial data sets for preclinical applications, it is easily incorporated into the existing DCS system, and it is flexible in that the source-detector separations and spatial scanning parameters are easily modifiable. Moreover, we have optimized the DCT experimental setup and image reconstruction algorithms for small animal leg imaging.
The paper is structured as follows: In Section 2, methods regarding the instrumentations, basic theory of DCS and DCT, comparison between a contact and a non-contact system, generation of simulated data sets from heterogeneous flow distribution, measurements of tissue phantoms made with a murine femur or a tube, and in vivo measurements on murine legs before and after allograft surgery are presented. In Section 3, we show the equivalent performance between a contact and a non-contact system, reconstruction results using simulated data, depth-resolved 3D DCT images of tissue phantoms, and in vivo temporal and spatial blood flow changes in murine allografts. A discussion of the results is presented in Section 4.
2.1 Diffuse correlation instrument
The diffuse correlation instrument generally consists of a light source, multiple photon-counting detectors and a correlator (Fig. 1(a) ). A continuous-wave, long-coherence-length (> 5m) diode laser working at 785 nm (DL785-120-SO, Crystal Laser, Reno, NV) is utilized as a light source. A contact probe containing optical fibers is placed in contact with the in vivo tissue of interest. Multi-mode fibers deliver near-infrared photons from the light source to the tissue. Photons propagated through tissue are collected from the tissue surface at different source-detector separations with single-mode fibers. Signals from each source-detector separation are detected by four photon-counting avalanche photodiodes (SPCM-AQ4C, Excelitas Technologies, Waltham, MA) and relayed to a 4-channel correlator board (Correlator.com, Bridgewater, NJ), which outputs the intensity temporal auto-correlation function.
2.2 Non-contact scanning module
To reconstruct three-dimensional blood flow images, diffuse correlation measurements at multiple source-detector combinations covering a large field of view are required. A non-contact scanning module that provides a large field of view, high positional accuracy, versatility of source-detector separation choices, and easy incorporation into an existing instrument is developed as an additional component added after the probe (Fig. 1(b)). Two identical lenses (LBF254-100-B, focal length f = 100 mm, Thorlabs, Newton, NJ) were utilized to image the probe onto the scanning plane. The lens placed near the probe collimated light from the multi-mode fiber and the second lens focused the collimated beam onto the scanning plane. A dual-axis scanning mirror galvanometer system (GVS012, Thorlabs, Newton, NJ) placed between the two lenses was used to redirect the collimated beam so that the probe could be imaged at different positions in the scanning plane. To minimize direct reflection from the scanning surface, crossed polarizers (APIR29-020, American Polarizers, Inc., Reading, PA) were placed before the source and detectors. Since the polarization of directly reflected light was set by the first polarizer, yet that of diffused light was randomized, the crossed polarizers greatly suppressed direct reflection while letting enough diffused signal pass. This scanning system had a capability to scan over an area of 30 mm x 30 mm with spatial resolution of 0.02 mm. The scanning system has no magnification (1:1), thus the source-detector separations were determined by those in the probe and could be easily modified for different applications. For this application, we utilized 3.0 mm and 4.2 mm source-detector separations.
2.3 Correlation diffusion equation
In dynamic light scattering, electric field auto-correlation function is defined as G 1(r,τ) ≡ ˂E*(r,t)·E(r,t + τ)˃, where E(r,t) is the electric field at position r and time t, τ is the delay time, * denotes complex conjugate and ˂ ˃ denotes integration over t. In highly scattering media (e.g., tissue) the propagation of G 1 can be described by the correlation diffusion equation, as follows ,
For geometries with simple boundary conditions such as infinite or semi-infinite geometries, the analytical solutions have been derived . The following is the analytic solution for the semi-infinite geometry:Eq. (2), is the unit vector pointing away from the turbid media boundary, ztr ≡ 1/ (μa + μ ś) is the depth where the source is considered isotropic, zb = 2(1 + Reff) / [3μ ś (1-Reff)] is the extrapolation distance from the boundary, Reff ≈-1.440n −2 + 0.710n −1 + 0.668 + 0.0636n is the effective reflection coefficient to account for the index mismatch between the tissue and air boundary and n is the refraction index of the tissue. ˂Δr 2˃ can be expressed into different forms based on the flow model used: for the case of random flow model, ˂Δr 2˃ = ˂V 2˃ τ 2, where ˂V 2˃ is the second moment of flow distribution. In the case of diffusive motion, such as Brownian motion, ˂Δr 2˃ = 6Dbτ, where Db is the Brownian diffusion coefficient. In practice, the Brownian motion model fits the observed auto-correlation decay of in vivo measurement better than random flow motion . Here, the Brownian motion model is assumed and αDb is referred to as the dynamic parameter of the medium that represents flow. For more complex geometries, numerical methods such as finite difference method (FDM) or finite element method (FEM) are usually utilized to solve the equation for G 1.
As is pointed out in the Section 2.1, the instrument provides measurements of the normalized intensity auto-correlation function g 2(τ); g 2(τ) ≡ ˂I(t)· I(t + τ)˃/˂I(t)˃2 where I(t) = |E(r,t)|2. The connection between the electric field auto-correlation and the intensity auto-correlation functions is provided by Siegert relation ,Eq. (3), β is a parameter primarily determined by the optical setup of the system, and is inversely proportional to the number of speckles in the detection area. For a contact system with a single mode fiber, β ≈ 0.5 as there are two independent linear polarization modes of speckle. For the proposed non-contact system, β ≈ 1.0 because the polarizers only allow one linearly polarized mode to pass.
2.4 Diffuse correlation tomography
With the diffusion correlation equation, the normalized temporal auto-correlation function g 1 can be calculated if the flow, αDb, is known. In practice, however, g 1 can be obtained from diffuse correlation measurement of g 2, while the flow distribution is unknown. The goal is to reconstruct the unknown distribution of flow based on the experimentally measured g 1 (i.e., solve the inverse problem).
To solve the inverse problem, each experimentally measured temporal auto-correlation function is divided into two parts, i.e., the contribution from homogeneous background and the perturbation from heterogeneity in flow, using a perturbation method. In this paper, the Rytov approximation was used to expand the temporal auto-correlation as g 1(r s,r d,τ) = g 1,0(r s,r d,τ)exp(ϕ(r s,r d,τ)), where g 1,0(r s,r d,τ) is the contribution from the homogeneous background and ϕ(r s,r d,τ) is the perturbation from heterogeneity of flow. r s is the position of the source and r d is the position of the detector. For computational simplicity, one representative data point from the correlation curve at a fixed delay time τ 0 is chosen for tomographic reconstruction, instead of calculating ϕ(r s,r d,τ) over the range of τs. Here, τ 0 which is the time when g 1,0, the calculated correlation curve of the background, decays to e −1 is chosen since it yields the optimal data set according to Zhou et al. .
By following the similar procedure in diffuse optical tomography  with the assumption that absorption and scattering coefficients are homogeneous and known, the perturbation, ϕ, can be expressed as follows:Eq. (4), i denotes the index for source-detector pair, and r is a position in the volume Ω. G 1,0(r si,r di,τ 0) is the solution to Eq. (1) with source at r si and detector at r di, and G 1,0(r si,r,τ 0) is the solution to Eq. (1) with source at r si and detector at r. H(r,r di,τ 0) is the Green’s function for Eq. (1) with source at r and detector at r di. To calculate the integral in Eq. (4), finite element mesh is used. First, the volume Ω is discretized into subspace: Ω1,…,Ωh,…, ΩL with nodes: r 1,…, r j,…, r N, where L is the total number of subspaces and N is the total number of nodes. In each subspace Ωh, the functions G, H and Δ(αDb) can be expressed as f h(r) = ∑j u j h(r,r j) f (r j), where f (r j) is function values at node r j and u j h(r,r j) is the shape function . Thus, Eq. (4) can be expressed in a matrix form as follows:Eq. (5), W ij establishes the link between the flow perturbation in the j th voxel, Δ(αDb(r j)), and the measurement perturbation at the i th source-detector pair. W ij is calculated as follows:Eq. (2). W is often called the sensitivity matrix or the weight matrix as it describes how sensitive the measurement is to the small changes in flow.
After weight matrix W is built, the last step to obtain the spatial distribution of flow change, Δ(αDb(r)), is to solve the linear equations listed in Eq. (5). Since the equations are usually ill-posed (i.e., the number of unknowns is much larger than that of linear independent equations), Tikhonov regularization is used to stabilize the matrix inversion:34].
2.5 Comparison between a non-contact and a contact system
To ensure that the quantitative performance of the non-contact system was equivalent to that of the contact system over the physiological flow range typically observed in vivo, both systems were directly compared following a method developed by Lin et al. . In this method, the Brownian motion of particles inside liquid phantom was increased by raising the temperature of the phantom gradually, which provided continuous flow changes over the physiological range. Concurrent measurements on a liquid phantom with μ a = 0.1 cm−1 and μ ś = 8.0 cm−1 were performed. The liquid phantom was made by mixing distilled water with nigrosin (Sigma-Aldrich, St. Louis, MO) as an absorbing agent and Intralipid (Fresenius Kabi, Uppsala, Sweden) as a scattering agent. The flow of the liquid phantom was changed gradually by increasing the temperature from 276.1K to 304.1K in 35 minutes. The temperature of the liquid phantom was increased by adding hot water to a tank surrounding the liquid phantom every 15 minutes and the temperature was monitored continuously using a thermometer (PS0614, Kent Scientific, Torrington, CT). The same integration time (2 s) and the source-detector separation (4.2 mm) were utilized for both the contact and the non-contact measurements as shown in Fig. 2 . Room light was turned off during the measurements.
The experimental data were fitted to the semi-infinite geometry solution in Eq. (2) to extract the flow, αDb. Regression analysis and Bland-Altman plot were utilized to compare αDb at different temperatures measured by the contact and the non-contact system.
2.6 DCT of simulated heterogeneous flow distributions
To test the performance of the DCT algorithm, image reconstructions were performed on a series of simulated data sets described in Fig. 3 and Table 1 . g 1 at multiple combinations of source-detector pairs from a three-dimensional map with heterogeneous distribution of αDb were simulated using a FEM-based solver of Eq. (1) modified from NIRFAST, a software package for diffuse optical tomography . The same mesh (x:-8.5 mm to 8.5 mm, y: −17.5 mm to 20.5 mm, z: 0 mm to 8.0 mm, total number of nodes: 136831, total number of elements: 673162) was used for all simulations. Figure 3(a) and 3(b) illustrates the side and the top view of the three-dimensional αDb map consisting of a cylinder of 1.2 mm diameter and 20 mm length simulating a long bone (Region 2) and the surrounding homogeneous medium simulating skeletal muscle (Region 1). Figure 3(c) shows the dense coverage of the source-detector pairs from scanning a source-detector pair of 3 mm separation with step size Δx of 0.35 mm and Δy of 0.70 mm. Noise was not added in the simulated data to test the inherent limitation of the linear reconstruction algorithm for the bone graft geometry.
In order to explore the effect of different μ ś between bone and muscle (20.0 cm−1 and 8.0 cm−1 ) in the reconstructed αDb images using the linear reconstruction algorithm based on the assumption of homogeneous background, four simulated data sets with various parameter assignments were produced as shown in Table 1. μ a was fixed as 0.1 cm−1 based on . Simulation #1 produced a data set similar to the experiment with ex vivo bone immersed in a liquid phantom described in Section 2.7, where αDb was lower in the bone (i.e., αDb = 0.25x10−8 cm2/s in region 2, αDb = 1.0x10−8 cm2/s in region 1), but μ ś of both regions were set to be 20.0 cm−1. Simulation #2 was the same as #1 except μ ś were set to be 8.0 cm−1 in both regions. Simulation #3 produced a data set with no αDb contrast, but with different μ ś in each region. Simulation #4 included both αDb contrast and μ ś differences.
For the DCT image reconstruction, a finite element mesh (x: −5.5 mm to 5.5 mm, y: −15.0 mm to 19.0 mm, z: 0 mm to 5.0 mm, total number of nodes: 80669, total number of elements: 374631) was used to reconstruct the weight matrix W. Note that this mesh was different from the mesh previously used to generate simulated data in terms of coverage and fineness of the mesh. For reconstruction of simulation #1, μ ś = 20.0 cm−1 was used in the reconstruction mesh. For reconstruction #2-4, μ ś = 8.0 cm−1 was assumed. The initial guess of αDb was 1.0x10−8 cm2/s for all cases.
After reconstructing a three-dimensional αDb map for each simulation, αDb values in region 2 were averaged to compare with the expected values (Table 1, αDb (2) column). Furthermore, relative blood flow (rBF) was calculated to quantify the changes before and after the αDb perturbation in a longitudinal study:
2.7 DCT of a tissue phantom with an ex vivo murine femur
To test the performance of the DCT reconstruction algorithm in the presence of realistic experimental noise, a tissue phantom was constructed by suspending an isolated murine femur in a liquid turbid medium. After harvesting a femur from a mouse, soft tissue debris attached to the femur was cleaned. Then the femur was suspended in a homogeneous liquid phantom with a thin wire going through the medullary cavity of the bone. The shaft of the femur has a cylindrical shape with about 1.2 mm diameter and 20 mm length. The dimension of liquid phantom was 100 mm x 70 mm x 60 mm in x, y, z direction. The longest axis of the femur was aligned in the y direction and the center of bone was 1.0 mm under the liquid phantom surface (i.e., z = 1.0 mm). The liquid phantom was made to match the absorption and scattering values of bone reported in literature , with μ a = 0.1 cm−1 and μ ś = 20.0 cm−1. This approach of matching optical properties between the bone and the surrounding medium simulated an idealized experimental situation for testing the linear reconstruction algorithm, which assumed homogeneous optical properties throughout the reconstruction volume.
The optical measurements were performed by scanning the surface of the turbid medium with the probe, covering the horizontal range x from −3.5 mm to 3.5 mm and the vertical range y from −12.6 mm to 15.6 mm. The step size for scanning was 0.35 mm in x direction and 0.70 mm in y direction. At each position, one auto-correlation curve was obtained with an integration time of 2 seconds. Data were taken at two source-detector separations, but only scanning data from 3.0 mm source-detector separation were utilized in the reconstruction for this preliminary study.
For the DCT image reconstruction, a finite element mesh was used (x: −4.5 mm to 4.5 mm, y: −13.5 mm to 16.5 mm, z: 0 mm to 6 mm, total number of nodes: 50060, total number of elements: 227788) to reconstruct the weight matrix W. For optical properties, μ a = 0.1 cm−1 and μ ś = 20.0 cm−1 were assumed. The background flow value was determined from a separate DCS measurement on the homogeneous liquid phantom before suspension of the femur, which was αDb = 0.82x10−8 cm2/s.
2.8 DCT of a tissue phantom with a tube
Although the use of an actual murine femur in the tissue phantom mimicked in vivo measurement geometry closely, it was difficult to modulate the flow within the bone. Therefore, another type of tissue phantom was constructed using a slightly turbid silicone rubber tube with an inner diameter of 1.0 mm (1972T4, McMaster Carr, Elmhurst, IL) instead of a femur. Two experiments were performed using the optical and flow parameters corresponding to simulation #3 and #4 (see Table 1). Lower flow inside the tube was achieved by adding 2% methyl cellulose to the liquid phantom to increase viscosity, thus decreasing Brownian motion of scatterers . Separate measurements on the liquid phantom yielded αDb = 0.90x10−8 cm2/s for liquid phantom without methyl cellulose and αDb = 0.19x10−8 cm2/s for phantom with 2% methyl cellulose.
The optical scan covered an area of x from −3.5 mm to 3.5 mm and y from −7.0 mm to 10.0 mm. The x scan step size was 0.35 mm and the y scan step size was 0.70 mm. One measurement was made at each position with integration time of 2 seconds.
For DCT reconstruction, a finite element mesh (x: −4.0 mm to 4.0 mm, y: −8 mm to 10.5 mm, z: 0 mm to 5 mm, total number of nodes: 52681, total number of elements: 244168) was used. The optical properties of the reconstruction regions were assumed to be homogeneous, with μ a = 0.1 cm−1 and μ ś = 8.0 cm−1.
2.9 DCT of in vivo murine legs before and after allograft surgery
For the in vivo experiment, three mice were scanned 5 days before and one week after an allograft surgery. All procedures in the animal study were approved by the University Committee on Animal Resource (UCAR) at the University of Rochester.
The details of the allograft surgery were described in Hoffman et al. [43–45 ]. To summarize, the allograft surgery was performed by removing a 5 mm diaphyseal segment of bone from the femur of the left hindlimb of a BALB/c mouse and replacing it with a decellularized bone segment from a C57BL/6 mouse (allograft). An intra-medullary pin through the host bones and the new allograft segment was placed to stabilize the graft and facilitate graft-host union.
Before DCT scanning, the mouse was placed on a warming pad and anesthetized with continuous isoflurane administration. Fur was removed from the scanning site 24 hours prior to baseline measurements using Nair (Church & Dwight Co., Inc., Princeton, NJ). Scanning field of view was adjusted such that the femur was located at the center of scanning area and the long axis was aligned with the y direction, as shown in Fig. 4 . For the second DCT scan, electric clippers were used to remove any regrown fur. A two dimensional scan was carried out covering an area of x from −3.5 mm to 3.5 mm and y from −7.0 mm to 10.0 mm. The x step size was 0.35 mm and the y step size was 0.70 mm. In total, measurements were made at 441 different positions with 2 source-detector separations and 2 seconds of integration time. After the measurement, the mouse was kept warm until it recovered from the anesthesia. The scanning time was approximately 18 minutes and the whole procedure took less than 40 minutes for each mouse. Only scanning data at 3.0 mm source-detector separation were used in the reconstruction.
For DCT image reconstruction, a finite element mesh was used (x: −4.0 mm to 4.0 mm, y: −8 mm to 10.5 mm, z: 0 mm to 5 mm, total number of nodes: 52681, total number of elements: 244168) to reconstruct the weight matrix, W. The optical properties of the reconstruction region were assumed to be homogeneous, with μ a = 0.1 cm−1 and μ ś = 8.0 cm−1. For each scan, first bulk αDb was fitted to the semi-infinite geometry solution in Eq. (2) at each source-detector pair, then bulk values were averaged over all pairs to yield a background homogeneous αDb needed for DCT reconstruction.
3.1 The performance of a non-contact system is equivalent to that of a contact-system
Typical normalized intensity auto-correlation functions measured by the non-contact system and contact system are shown in Fig. 5(a) . Note that experimental β of the non-contact and the contact system is close to 1.0 and 0.5 respectively, as mentioned in Section 2.3. A linear relationship between αDb measured by the contact system and non-contact system was observed as shown in Fig. 5(b). Linear regression yielded a slope of 0.99 with high coefficient of determination R 2 = 0.87 (p < 0.0001). The agreement between the two methods can be assessed by a Bland-Altman plot, constructed by plotting the average of the two variables along the x-axis and the difference between the two variables along the y-axis , as shown in Fig. 5(c). 95% limits of agreement were computed by the mean difference 1.96 standard deviation of the difference (dotted line). We can conclude that the contact and the non-contact systems generally agree since 96% of data lies within 95% limits of agreement.
From these results, it can be concluded that the performance of the non-contact system is comparable to the contact system. Lin et al.  also reported similar results comparing contact and non-contact systems for much larger source-detector separation (26 mm) than used in this study (4.2 mm).
3.2 Three-dimensional flow image reconstruction of simulated data sets
In Fig. 6(a) , three-dimensional images are shown as a series of 8 mm x 30 mm image slices at different depths (z) of 0.5 mm interval. The top row images are the expected flow distributions. The middle and the bottom row images are the reconstructed αDb images from simulation #1 and #2 respectively. The reconstructed images show that the extent of flow heterogeneity is larger than the expected size (i.e., 1.2 mm diameter) and the flow value inside the cylinder is underestimated (0.91x10−8 cm2/s vs expected 0.25x10−8 cm2/s). This is in line with the observations from Zhou et al.  and other diffuse optical tomography results  that the point spread function of the diffuse optical techniques is large and the reconstructed values are usually underestimated. Interestingly, the depth of the cylinder in simulation #1 with global μ ś = 20.0 cm−1 was underestimated (0.57 mm vs. expected 1.0 mm), yet that of simulation #2 with global μ ś = 8.0 cm−1 was overestimated (1.24 mm). The reconstructed images from simulation #3 and #4 (not shown) looked similar to the z-distribution of simulation #2.In Fig. 6(b), representative αDb slices at z = 1.0 mm from four simulations are shown side by side to demonstrate the effect of μ ś differences between regions. When no difference was present, reconstructed αDb in the cylindrical region was lower than the background as expected (simulation #2, mean αDb in region 2 = 0.88x10−8 cm2/s). However, μ ś differences in the absence of αDb contrast between two regions (simulation #3) resulted in the apparent reconstructed αDb increase in the cylindrical region. When both μ ś differences and αDb contrast were present (simulation #4), reconstructed αDb was lower in the cylindrical region, but the contrast (0.94x10−8 cm2/s) was less than that of simulation #2. For the longitudinal study, the accurate quantification of flow changes from the baseline may still be possible even though the accuracy of absolute αDb is compromised due to systematic error of the reconstruction method. Therefore, the distribution of relative flow was computed by dividing reconstructed αDb map of simulation #4 with that of simulation #3, and compared with the expected rBF map in Fig. 6(c). rBF in the cylindrical region was 0.78 as opposed to the expected rBF of 0.25.
3.3 Three-dimensional flow image reconstruction of a tissue phantom with a murine femur
Figure 7 shows the reconstructed αDb images in different layers of the tissue phantom at 0.5 mm increments from the surface of the liquid phantom. The photograph on the left shows the position of the bone with respect to the scanning region. At z = 0.5 and 1.0 mm, αDb distribution with a lower value than the background is observed at the approximate region that corresponds to the spatial extent of the bone in x-y plane, with the narrower width at the center and wider width at the end of the bone. As the bone is solid and, therefore, has a much smaller diffusion coefficient than the liquid phantom, it is expected to exhibit lower flow than the surrounding liquid phantom. The average reconstructed z depth of the bone region was 0.7 mm, which was smaller than the expected depth of 1.0 mm.
3.4 Three-dimensional flow image reconstructions of a tissue phantom with a tube
Figure 8 shows the experimental analogue of Fig. 6(c). rBF images are presented in different layers of the tissue phantom at 0.5 mm increments from the surface of the liquid phantom inthe presence of μ ś differences between regions. Lower relative flow (rBF = 0.68) in the tube region was observed, but the extent was underestimated when compared to the expected relative flow (rBF = 0.22).
3.5 Three-dimensional blood flow image reconstructions of in vivo murine leg before and after allograft surgery
In Fig. 9 , rBF images of all three mice at different depths (0 mm, 0.5 mm, 1 mm, 1.5 mm, 2.0 mm and 2.5 mm below the scanning surface) are shown. rBF images were calculated by dividing reconstructed αDb from measurements taken one week after surgery with that from measurements before surgery. The approximate location of the femur and the graft is outlined with a black dashed line and a solid red line respectively. For every mouse, overall elevation of blood flow was observed compared to the pre-surgical state (1.5 ~2 times). In addition, localized blood flow changes near the graft were clearly captured for all three mice imaged, with the maximum rBF value ranging from 3 to 5.
In this study, we demonstrated a galvanometer-based non-contact diffuse correlation tomography system that can provide large spatial data sets with easy instrumental implementation and flexibility, which in turn provides more spatially resolved information in three-dimensional blood flow images of murine femurs. Previous DCT works have been performed with a limited number of spatial data sets for rat brain imaging [33, 34 ]. Recently, new techniques based on speckle-contrast have been emerging as an alternative for deep-tissue blood flow imaging [48–51 ]. This can provide large spatial data sets using a CCD-based system, which can yield resolution equivalent to the traditional diffuse correlation tomographic approaches. However, these works were focused on applications in the human scale and the full capability of improving resolution has not been explicitly demonstrated as of yet.
Compared to laser Doppler imaging or speckle contrast imaging, the large data set DCT approach is ideally suited for monitoring the blood flow of the murine bone graft model because of the requirements for relatively deep penetration depth and a large field of view. The murine femur is approximately 1-2 mm deep with respect to the tissue surface, and the entire hindlimb is approximately 5 mm thick. The location of the femur requires relatively deep tissue penetration of light. The length of the femur is approximately 15 - 20 mm with the implanted graft being 5 mm long. This configuration requires not only depth penetration, but also a relatively large field of view to cover the whole femur region. We can provide this large field of view through the galvanometer-based scanning system.
At the present time, the utilized DCT instrumentation and algorithm are not fully optimized and need further improvements. The optical throughput of the current scanning module was lower than that of the contact system due to low transmission of polarizers (30% at 785nm) and insufficient optimization of optical components for light collection efficiency. While the light loss due to crossed polarizers is unavoidable in our scanning design (i.e., 9% transmission after two polarizers), the light throughput in other optical components can be improved with a better design. In the experiments presented here, no optical components were added to equalize the signal-to-noise at different source-detector separations, resulting in superior signal-to-noise ratio for 3 mm separation compared to that for 4.2 mm. For this preliminary study, we focused on using only 3 mm separation data to simplify the analysis. Testing new optical designs to improve the light throughput and to equalize the signal-to-noise at multiple separations is ongoing.
The most notable discrepancy between reconstructed and expected flow distribution using simulated data was the underestimation of flow contrast between the heterogeneity and the surrounding medium. Our current algorithm was based on a linear Rytov approximation which usually underestimates flow values . Based on simulations (not shown), inclusion of more source-detector separations in the DCT reconstruction can improve the accuracy of reconstructed flow values. Another discrepancy was the depth position of the heterogeneity which may be due to the use of constant Tikhonov regularization. This can be improved with spatially variant regularization . In the future, the strategy to optimally combine data from multiple separations with varying signal-to-noise ratios will be developed by inspecting sensitivity matrices of various optical and flow parameters, in addition to incorporating an iterative reconstruction scheme.
Another improvement for DCT reconstruction to be made is the incorporation of heterogeneous distribution of absorption and reduced scattering within the algorithm. In the current reconstruction, absorption and scattering were assumed to be homogeneous throughout. In the ex vivo murine femur tissue phantom experiment, absorption and scattering of the liquid phantom was matched with reported values of bone in literature to create an ideal experimental condition. For the animal study, absorption and scattering coefficients were assumed to match those of muscle, and the potential scattering coefficient differences between the muscle and the bone was not incorporated in the reconstruction algorithm. The impact of this assumption on the quantification of absolute and relative blood flow was studied using a series of numerical simulations and tissue phantom experiments implementing scattering coefficient differences between the muscle and the bone (see Sections 2.6 and 3.2). We showed that relative flow was more robust and reliable compared to the absolute flow, which was susceptible to the optical property differences between the bone and the surrounding tissue. One way to improve this limitation in absolute flow recovery is to use mouse-specific anatomical information from other imaging modalities such as micro-CT, and to assign specific absorption and scattering coefficients from literature in different tissues such as muscle and bone. Another way is to combine DCT with diffuse optical tomography (DOT) and assign the three-dimensional absorption and scattering coefficient map based on DOT reconstruction before performing DCT reconstruction.
Albeit current limitations, we demonstrated superior lateral resolution in a realistic scale tissue phantom with an immersed murine femur using a large spatial data set. In addition, DCT images from before and after allograft surgery revealed overall elevation of blood flow level after surgery and spatially localized flow increases within the femur and/or the surrounding tissue. The ability to reliably and longitudinally monitor blood flow in the murine segmental bone defect model is important as the therapeutic efficacy of various tissue engineering approaches to enhance allograft healing is being assessed using this model [43–45, 53 ]. In particular, the number and size of blood vessels during bone repair and regeneration have been extensively quantified using end-point assays such as histology or micro-CT with a lead-chromate contrast agent . However, blood flow, which dictates the rate of delivery of oxygen, nutrients, growth factors and circulating cells necessary for healing, is not readily available for quantification. Despite the difficulties, blood flow in the region surrounding a bone injury is thought to be increased due to inflammation and metabolic demands in the early phase of healing . In the later phase of healing, blood flow is expected to increase within the bone graft after the establishment of new penetrating blood vessels from existing vessels . Our preliminary DCT images taken one week after allograft surgery suggest that blood flow needed for healing may be increased in the tissue surrounding the bone graft, considering the limited depth resolution of current analysis. The longitudinal monitoring showing spatial evolution of blood flow in the same mouse at different healing phases may reveal the interplay between the surrounding tissue and graft, and may yield quantitative parameters to predict therapeutic efficacy in the future.
In summary, we have constructed a scanning non-contact DCT system with galvanometer mirrors and verified that our non-contact system performed equivalently to a contact system. We also demonstrated the feasibility of three-dimensional flow imaging with our new DCT system from a tissue phantom with a femur and a tube, and in vivo murine legs. Overall blood flow after 1 week post-allograft surgery on murine legs increased by 1.5 to 2 times from the baseline (i.e., before surgery) blood flow in all three mice. In addition, higher relative blood flow (maximum 3 to 5 times from the baseline) was observed in a focused spatial region near the graft in each mouse. Our system demonstrated a potential for longitudinal study of blood flow changes during bone graft healing.
We thank Emmanuel Mannoh for help with in vivo measurements, Tanja Dragojević, Martina Giovannella, Marco Pagliazzi, Clara Gregori Pla and Hari M. Varma for the helpful discussion in devising phantom experiments, Jingxuan Ren, Gabriel Ramirez for their support. We also thank Michael Thullen and Dr. Hani Awad for musculoskeletal research support. This research was supported by NIH K99/R00-CA126187 (Choe), NIH T32 AR053459 (Hoffman), NIH R01 AR064200 (Benoit), NYSTEM IDEA, C028099 (Benoit), Fundació Cellex Barcelona (Durduran), Ministerio de Economía y Competitividad (PHOTOSTROKE) (Durduran), l'Obra Social La Caixa (MedLlumBCN) (Durduran) and LASERLAB Europe III (Bioptichal) (Durduran).
References and Links
1. A. S. Greenwald, S. D. Boden, V. M. Goldberg, Y. Khan, C. T. Laurencin, R. N. Rosier, and American Academy of Orthopaedic Surgeons. The Committee on Biological Implants, “Bone-graft substitutes: facts, fictions, and applications,” J. Bone Joint Surg. Am. 83-A(Suppl 2 Pt 2), 98–103 (2001). [PubMed]
3. F. J. Hornicek, M. C. Gebhardt, W. W. Tomford, J. I. Sorger, M. Zavatta, J. P. Menzner, and H. J. Mankin, “Factors affecting nonunion of the allograft-host junction,” Clin. Orthop. Relat. Res. 382, 87–98 (2001). [CrossRef] [PubMed]
6. R. E. Tomlinson and M. J. Silva, “Skeletal blood flow in bone repair and maintenance,” Bone Res. 1(4), 311–322 (2013). [CrossRef]
8. M. E. Phelps, J. C. Mazziotta, H. R. Schelbert, R. A. Hawkins, and J. Engel, “Clinical PET-what are the issues,” J. Nucl. Med. 26, 1353–1358 (1985).
9. P. Lindholm, E. Sutinen, V. Oikonen, K. Mattila, M. Tarkkanen, M. Kallajoki, H. Aro, T. Böhling, A. Kivioja, I. Elomaa, and H. Minn, “PET imaging of blood flow and glucose metabolism in localized musculoskeletal tumors of the extremities,” Nucl. Med. Biol. 38(2), 295–300 (2011). [CrossRef] [PubMed]
11. N. Lassau, C. Paturel-Asselin, J. M. Guinebretiere, J. Leclère, S. Koscielny, A. Roche, S. Chouaib, and P. Peronneau, “New hemodynamic approach to angiogenesis: color and pulsed Doppler ultrasonography,” Invest. Radiol. 34(3), 194–198 (1999). [CrossRef] [PubMed]
12. T. S. Padayachee, F. J. Kirkham, R. R. Lewis, J. Gillard, M. C. Hutchinson, and R. G. Gosling, “Transcranial measurement of blood velocities in the basal cerebral arteries using pulsed Doppler ultrasound: a method of assessing the Circle of Willis,” Ultrasound Med. Biol. 12(1), 5–14 (1986). [CrossRef] [PubMed]
13. N. M. Tole, Basic Physics of Ultrasonographic Imaging (World Health Organization Press, Switzerland, 2005).
14. T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73(7), 076701 (2010). [CrossRef]
15. C. Cheung, J. P. Culver, K. Takahashi, J. H. Greenberg, and A. G. Yodh, “In vivo cerebrovascular measurement combining diffuse near-infrared absorption and correlation spectroscopies,” Phys. Med. Biol. 46(8), 2053–2065 (2001). [CrossRef] [PubMed]
16. G. Yu, T. Durduran, G. Lech, C. Zhou, B. Chance, E. R. Mohler 3rd, and A. G. Yodh, “Time-dependent blood flow and oxygenation in human skeletal muscles measured with noninvasive near-infrared diffuse optical spectroscopies,” J. Biomed. Opt. 10(2), 024027 (2005). [CrossRef] [PubMed]
17. G. Yu, T. Durduran, C. Zhou, H. W. Wang, M. E. Putt, H. M. Saunders, C. M. Sehgal, E. Glatstein, A. G. Yodh, and T. M. Busch, “Noninvasive monitoring of murine tumor blood flow during and after photodynamic therapy provides early assessment of therapeutic efficacy,” Clin. Cancer Res. 11(9), 3543–3552 (2005). [CrossRef] [PubMed]
18. U. Sunar, H. Quon, T. Durduran, J. Zhang, J. Du, C. Zhou, G. Yu, R. Choe, A. Kilger, R. Lustig, L. Loevner, S. Nioka, B. Chance, and A. G. Yodh, “Noninvasive diffuse optical measurement of blood flow and blood oxygenation for monitoring radiation therapy in patients with head and neck tumors: a pilot study,” J. Biomed. Opt. 11(6), 064021 (2006). [CrossRef] [PubMed]
19. C. Zhou, R. Choe, N. Shah, T. Durduran, G. Yu, A. Durkin, D. Hsiang, R. Mehta, J. Butler, A. Cerussi, B. J. Tromberg, and A. G. Yodh, “Diffuse optical monitoring of blood flow and oxygenation in human breast cancer during early stages of neoadjuvant chemotherapy,” J. Biomed. Opt. 12(5), 051903 (2007). [CrossRef] [PubMed]
20. C. Zhou, S. A. Eucker, T. Durduran, G. Yu, J. Ralston, S. H. Friess, R. N. Ichord, S. S. Margulies, and A. G. Yodh, “Diffuse optical monitoring of hemodynamic changes in piglet brain with closed head injury,” J. Biomed. Opt. 14(3), 034015 (2009). [CrossRef] [PubMed]
22. T. Durduran, “Non-invasive measurements of tissue hemodynamics with hybrid diffuse optical methods,” (University of Pennsylvania, 2004).
23. R. C. Mesquita, N. Skuli, M. N. Kim, J. Liang, S. Schenkel, A. J. Majmundar, M. C. Simon, and A. G. Yodh, “Hemodynamic and metabolic diffuse optical monitoring in a mouse model of hindlimb ischemia,” Biomed. Opt. Express 1(4), 1173–1187 (2010). [CrossRef] [PubMed]
24. C. Menon, G. M. Polin, I. Prabakaran, A. Hsi, C. Cheung, J. P. Culver, J. F. Pingpank, C. S. Sehgal, A. G. Yodh, D. G. Buerk, and D. L. Fraker, “An integrated approach to measuring tumor oxygen status using human melanoma xenografts as a model,” Cancer Res. 63(21), 7232–7240 (2003). [PubMed]
25. E. M. Buckley, N. M. Cook, T. Durduran, M. N. Kim, C. Zhou, R. Choe, G. Yu, S. Schultz, C. M. Sehgal, D. J. Licht, P. H. Arger, M. E. Putt, H. H. Hurt, and A. G. Yodh, “Cerebral hemodynamics in preterm infants during positional intervention measured with diffuse correlation spectroscopy and transcranial Doppler ultrasound,” Opt. Express 17(15), 12571–12581 (2009). [CrossRef] [PubMed]
26. N. Roche-Labarbe, S. A. Carp, A. Surova, M. Patel, D. A. Boas, P. E. Grant, and M. A. Franceschini, “Noninvasive optical measures of CBV, StO(2), CBF index, and rCMRO(2) in human premature neonates’ brains in the first six weeks of life,” Hum. Brain Mapp. 31(3), 341–352 (2010). [CrossRef] [PubMed]
27. P. Zirak, R. Delgado-Mederos, J. Martí-Fàbregas, and T. Durduran, “Effects of acetazolamide on the micro- and macro-vascular cerebral hemodynamics: a diffuse optical and transcranial doppler ultrasound study,” Biomed. Opt. Express 1(5), 1443–1459 (2010). [CrossRef] [PubMed]
28. T. Durduran, G. Yu, M. G. Burnett, J. A. Detre, J. H. Greenberg, J. Wang, C. Zhou, and A. G. Yodh, “Diffuse optical measurement of blood flow, blood oxygenation, and metabolism in a human brain during sensorimotor cortex activation,” Opt. Lett. 29(15), 1766–1768 (2004). [CrossRef] [PubMed]
29. G. Yu, T. F. Floyd, T. Durduran, C. Zhou, J. Wang, J. A. Detre, and A. G. Yodh, “Validation of diffuse correlation spectroscopy for muscle blood flow with concurrent arterial spin labeled perfusion MRI,” Opt. Express 15(3), 1064–1075 (2007). [CrossRef] [PubMed]
30. T. Durduran, C. Zhou, E. M. Buckley, M. N. Kim, G. Yu, R. Choe, J. W. Gaynor, T. L. Spray, S. M. Durning, S. E. Mason, L. M. Montenegro, S. C. Nicolson, R. A. Zimmerman, M. E. Putt, J. Wang, J. H. Greenberg, J. A. Detre, A. G. Yodh, and D. J. Licht, “Optical measurement of cerebral hemodynamics and oxygen metabolism in neonates with congenital heart defects,” J. Biomed. Opt. 15(3), 037004 (2010). [CrossRef] [PubMed]
31. S. A. Carp, G. P. Dai, D. A. Boas, M. A. Franceschini, and Y. R. Kim, “Validation of diffuse correlation spectroscopy measurements of rodent cerebral blood flow with simultaneous arterial spin labeling MRI; towards MRI-optical continuous cerebral metabolic monitoring,” Biomed. Opt. Express 1(2), 553–565 (2010). [CrossRef] [PubMed]
32. M. N. Kim, T. Durduran, S. Frangos, B. L. Edlow, E. M. Buckley, H. E. Moss, C. Zhou, G. Yu, R. Choe, E. Maloney-Wilensky, R. L. Wolf, M. S. Grady, J. H. Greenberg, J. M. Levine, A. G. Yodh, J. A. Detre, and W. A. Kofke, “Noninvasive measurement of cerebral blood flow and blood oxygenation using near-infrared and diffuse correlation spectroscopies in critically brain-injured adults,” Neurocrit. Care 12(2), 173–180 (2010). [CrossRef] [PubMed]
33. J. P. Culver, T. Durduran, D. Furuya, C. Cheung, J. H. Greenberg, and A. G. Yodh, “Diffuse optical tomography of cerebral blood flow, oxygenation, and metabolism in rat during focal ischemia,” J. Cereb. Blood Flow Metab. 23(8), 911–924 (2003). [CrossRef] [PubMed]
34. C. Zhou, G. Yu, D. Furuya, J. Greenberg, A. Yodh, and T. Durduran, “Diffuse optical correlation tomography of cerebral blood flow during cortical spreading depression in rat brain,” Opt. Express 14(3), 1125–1144 (2006). [CrossRef] [PubMed]
36. Y. Lin, C. Huang, D. Irwin, L. He, Y. Shang, and G. Yu, “Three-dimensional flow contrast imaging of deep tissue using noncontact diffuse correlation tomography,” Appl. Phys. Lett. 104(12), 121103 (2014). [CrossRef] [PubMed]
38. P. A. Lemieux and D. J. Durian, “Investigating non-Gaussian scattering processes by using nth-order intensity correlation functions,” J. Opt. Soc. Am. A 16(7), 1651–1664 (1999). [CrossRef]
39. M. A. O'Leary, “Imaging with diffuse photon density waves,” (University of Pennsylvania, 1996).
40. T. R. Chandrupatla and A. D. Belegundu, Introduction to finite elements in engineering (Prentice Hall, Upper Saddle River, N.J., 2002).
41. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. 25(6), 711–732 (2009). [CrossRef] [PubMed]
43. M. D. Hoffman, C. Xie, X. Zhang, and D. S. Benoit, “The effect of mesenchymal stem cells delivered via hydrogel-based tissue engineered periosteum on bone allograft healing,” Biomaterials 34(35), 8887–8898 (2013). [CrossRef] [PubMed]
44. M. D. Hoffman, A. H. Van Hove, and D. S. W. Benoit, “Degradable hydrogels for spatiotemporal control of mesenchymal stem cells localized at decellularized bone allografts,” Acta Biomater. 10(8), 3431–3441 (2014). [CrossRef] [PubMed]
45. M. D. Hoffman and D. S. W. Benoit, “Emulating native periosteum cell population and subsequent paracrine factor production to promote tissue engineered periosteum-mediated allograft healing,” Biomaterials 52, 426–440 (2015). [CrossRef] [PubMed]
46. M. Bland, An Introduction to Medical Statistics (Oxford University Press, Oxford; New York, 2000).
47. J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, and A. G. Yodh, “Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,” Med. Phys. 30(2), 235–247 (2003). [CrossRef] [PubMed]
48. H. M. Varma, C. P. Valdes, A. K. Kristoffersen, J. P. Culver, and T. Durduran, “Speckle contrast optical tomography: A new method for deep tissue three-dimensional tomography of blood flow,” Biomed. Opt. Express 5(4), 1275–1289 (2014). [CrossRef] [PubMed]
49. C. P. Valdes, H. M. Varma, A. K. Kristoffersen, T. Dragojevic, J. P. Culver, and T. Durduran, “Speckle contrast optical spectroscopy, a non-invasive, diffuse optical method for measuring microvascular blood flow in tissue,” Biomed. Opt. Express 5(8), 2769–2784 (2014). [CrossRef] [PubMed]
50. R. Bi, J. Dong, and K. Lee, “Multi-channel deep tissue flowmetry based on temporal diffuse speckle contrast analysis,” Opt. Express 21(19), 22854–22861 (2013). [PubMed]
52. B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Österberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38(13), 2950–2961 (1999). [CrossRef] [PubMed]
53. J. A. Inzana, D. Olvera, S. M. Fuller, J. P. Kelly, O. A. Graeve, E. M. Schwarz, S. L. Kates, and H. A. Awad, “3D printing of composite calcium phosphate and collagen scaffolds for bone regeneration,” Biomaterials 35(13), 4026–4034 (2014). [CrossRef] [PubMed]
54. C. L. Duvall, W. R. Taylor, D. Weiss, and R. E. Guldberg, “Quantitative microcomputed tomography analysis of collateral vessel development after ischemic injury,” Am. J. Physiol. Heart Circ. Physiol. 287(1), H302–H310 (2004). [CrossRef] [PubMed]