Non-invasive diffuse correlation tomography reveals spatial and temporal blood flow differences in murine bone grafting approaches.

Longitudinal blood flow during murine bone graft healing was monitored non-invasively using diffuse correlation tomography. The system utilized spatially dense data from a scanning set-up, non-linear reconstruction, and micro-CT anatomical information. Weekly in vivo measurements were performed. Blood flow changes in autografts, which heal successfully, were localized to graft regions and consistent across mice. Poor healing allografts showed heterogeneous blood flow elevation and high inter-subject variabilities. Allografts with tissue-engineered periosteum showed responses intermediate to both autografts and allografts, consistent with healing observed. These findings suggest that spatiotemporal blood flow changes can be utilized to differentiate the degree of bone graft healing.


Introduction
Critical-sized bone defects cannot heal without intervention [1]. Annually, more than 500,000 U.S. and 2.2 million global orthopedic procedures are carried out to treat critical-sized bone defects [2,3]. Among those procedures, autograft transplantation utilizes healthy endogenous bone and soft tissue from a non-load bearing region of the skeleton. Autografts often achieve complete healing; however, the procedure is limited by tissue availability and associated donor site morbidity [2,3]. Alternatively, allograft transplantation, the "gold-standard" treatment for critical-sized bone defects [2,4,5], utilizes processed cadaveric bone to repair bone defects. To reduce immunogenicity, allografts are completely devitalized, during which the periosteum, a thin tissue layer covering the outer bone surface, which is important for complete bone healing [6], is lost. As a result, allografts exhibit poor host osteointegration with a 60%, 10-year post-implantation failure rate [7][8][9]. To improve allograft integration and healing, tissue engineered (T.E.) allografts have been developed and studied. The T.E. approaches deliver periosteum mimicking cells and/or growth factors to the graft site, which include the use of bone morphogenetic proteins, mesenchymal stem cells (MSCs), and/or osteoprogenitor cells [10][11][12][13]. To assess the effectiveness of allografts augmented with tissue engineered periostea, a murine segmental femoral defect model (graft model) is often used in pre-clinical studies. In the graft model, the healing outcome and the physiological changes during healing are assessed, such as graft vascularization, bone callus formation and bone strength [5].
Vascularization is recognized as an especially critical process in bone healing, as blood supplies the necessary nutrients, circulating cells, growth factors, and oxygen to the graft site [14]. Overall blood supply to grafts depends both on intravascular volume and blood flow. Blood flow can change rapidly in existing vasculature of native bones and tissues surrounding grafts due to inflammation and metabolic demands in the early phase of healing [14]. Characterization of longitudinal blood flow changes, as an important step to study graft vascularization, has the potential not only to advance our understanding of graft healing, but to establish a standard to compare the healing process in different graft types. However, the measurement of blood flow in bone is often confounded by high bone density as well as the heterogeneity of the intraosseous blood supply [15]. A few methods to measure bone blood flow are available, such as radioactive or fluorescent microsphere techniques [16,17], magnetic resonance imaging (MRI), positron emission tomography (PET) [15], and laser Doppler flowmetry [18,19]. Nevertheless, the application of these techniques for longitudinal studies of blood flow in bone has proven difficult due to invasiveness (microspheres), limited penetration depth (laser Doppler) or cost (MRI and PET) [15,19,20]. As a result, systematic comparison of spatial and temporal blood flow changes in different graft types has not yet been performed.
Diffuse correlation spectroscopy (DCS) and tomography (DCT) are deep-tissue, noninvasive techniques that can quantify blood flow index, which is sensitive to microvasculature [21]. DCS and DCT are developed by expanding the theory of dynamic light scattering, which mainly accounts for a single scattering event, to the domain where multiple scattering dominates, which is the case after near-infrared light (650 nm -950 nm) propagates more than a few millimeters within highly scattering tissue [22,23]. With DCS/DCT, the fluctuation of collected light signal, quantified by the intensity auto-correlation function, is linked to the motion of scatterers (i.e., red blood cells) [21].The measurements are noninvasive, inexpensive and easy to perform. The measurement method to infer average flow index from the diffuse light intensity auto-correlation function is referred to as DCS and the 3-dimensional (3D) imaging method is called DCT.
DCS has been widely applied to monitor the blood flow in tumors, brain and skeletal muscle in animals and humans in vivo [24][25][26][27][28][29][30] [29,42]. In most cases, the relative change of blood flow with respect to the baseline has been validated. However, careful calibration studies can be performed to convert the blood flow index to a quantitative measure of perfusion [43]. DCT, which provides 3D relative blood flow images, has been applied to monitor 3D blood flow in rat brain [44,45] and human breasts [46]. Previously, we reported a non-contact scanning DCT system designed for murine graft study and showed its capability to capture localized flow changes between before and 1 week after graft transplantation in mice with allografts [47].
In this study, we report improvement on the DCT system as well as the results of longitudinal blood flow changes in mice with three different graft types; i.e., autograft, allograft and T.E. allograft with MSCs. The paper is structured as follows: In Section 2, methods regarding the instrumentation, a nonlinear DCT reconstruction, simulation set-up to compare a linear and a non-linear DCT reconstruction method, in vivo mouse measurement protocol and generation of mesh from micro-CT images are presented. In Section 3, we show the performance comparison between linear and non-linear DCT reconstruction methods, a representative example of the longitudinal 3D blood flow changes in a mouse with an autograft, comparison of blood flow changes in three representative mice from different groups and the averaged blood flow changes in the graft region. Finally, a discussion of the results is presented in Section 4.

Instrumentation
A typical diffuse correlation instrument consists of a light source, photon-counting avalanche photodiodes (APDs) and a correlator, as shown on the left side of Fig. 1. The long-coherencelength laser source at 785 nm is delivered to the tissue using a multi-mode fiber. Multiplescattered light signals are collected using single-mode fibers at several millimeters away from the source. Photon-counting avalanche photodiodes (SPCM-AQ4C, Excelitas Technologies, Waltham, MA) are utilized to detect the light intensity and relay the data to a correlator (Flex030EM, Correlator.com, Bridgewater, NJ), which outputs the intensity temporal autocorrelation function g 2 .
To collect spatially dense data for DCT, a non-contact scanning module was added to the typical DCS instrument as shown on the right of Fig. 1. A fiber collimator is employed to collimate the beam coming out of the source fiber. Then, the collimated beam goes through one arm of a polarizing beamsplitter (CM1-PBS252, Thorlabs, Newton, NJ), is redirected by two galvanometer mirrors and projected on the tissue surface using a 50 mm doublet lens (AC254-050-B, Thorlabs, Newton, NJ). The typical beam spot size and the light power on the tissue is 0.3 mm and less than 1 mW, respectively. For the data collection, the detection fibers are imaged onto the tissue surface at known separations from the source through the other arm of the polarizing beamsplitter using two identical lenses. With this design, the sourcedetector separations are maintained while the galvanometer mirrors scan the source-detector pattern in two dimensions. The current design is improved upon our previous system [47], and it has better light throughput. The scanning module is similar to the scanning module design of the laminar optical tomography (LOT) system [48]; however, this particular noncontact scanning module enables measurements at much larger source-detector separations and has a bigger field of view than the LOT system.

Non-linear image reconstruction
Using the diffuse correlation instrument, the normalized intensity auto-correlation function g 2 (r,τ) ≡ ˂I(r,t)·I(r,t + τ)˃/˂I(r,t)˃ 2 is measured, where τ is the delay time, I(r,t) in the measured intensity over time t at position r, and I(r,t) = |E(r,t)| 2 with E(r,t) representing the electric field vector. The fundamental quantity under investigation in optical coherence theory is the electric field temporal auto-correlation function, G 1 (r,τ) ≡ ˂E*(r,t)·E(r,t + τ)˃, where * denotes complex conjugate and < > denotes integration over time. The Siegert relation provides the connection between the normalized intensity auto-correlation and the normalized electric field auto-correlation functions [49], where the normalized electric field temporal auto-correlation is g 1 (τ) ≡ G 1 (r,τ)/ G 1 (r,0). In Eq.
(1), β is a parameter primarily determined by the optical set-up of the system, and is inversely proportional to the number of speckles in the detection area [49].
In highly scattering media (e.g., tissue), the correlation diffusion equation provides the theoretical model to describe the propagation of electric field temporal auto-correlation, G 1 (r,τ), as follows: where μ a and μ ś are the absorption and reduced scattering coefficient respectively, ν is the speed of light in the medium, D ≡ ν/(3(μ a + μ ś) ) is the photon diffusion coefficient, κ 0 = 2π/λ is the wave vector of the photon where λ is the wavelength in the tissue, α is the fraction of moving scatterers and ˂Δr 2 (r,τ)˃ is the mean-squared displacement of scatterers, Sδ 3 (r-r s ) is a point source located at r s with the source strength S. In Eq. (2), flow is characterized using ˂Δr 2 (r,τ)˃. To simplify this expression, a Brownian motion model is usually utilized [21], which gives ˂Δr 2 ˃ = 6D b τ and D b is the Brownian diffusion coefficient. As a result, the parameter, αD b , is referred to as the blood flow index.
To solve for G 1 (r,τ) in Eq. (2), the following partial-flux boundary condition is used [21, 50]: where n is the unit vector pointing away from the scattering media boundary and 0636n is the effective reflection coefficient to account for the index mismatch between the tissue and air boundary, and n is the refraction index ratio of tissue and air. Analytical solutions to the correlation diffusion equation in a simple geometry, for example infinite, semi-infinite, or layered geometry, have been derived and widely used in various studies [21]. In more general cases, the correlation diffusion equation needs to be solved numerically using finite difference method (FDM) or finite element method (FEM).
With the experimental measurements, g 2 is obtained for multiple source-detector pairs and the corresponding g 1,m can be calculated using Eq. (1). DCT utilizes the experimentally obtained g 1,m to reconstruct the flow distribution, αD b (r): this process is called solving the inverse problem. The inverse problem is constructed to minimize a cost function of the Rytov fashion [51], as follows: In Eq. (4), N sd is the total number of source-detector pairs utilized in the measurements, i represents one of the source-detector pairs, r si is the corresponding source position, r di is the detector position, g 1,m (r si ,r di ,τ) is the measured auto-correlation function with source located at r si and detector at r di , and g 1,c (r si ,r di ,τ) is the calculated auto-correlation function based on an estimate of flow distribution using Eq. (2). As the auto-correlation function spans over a wide range of τ values, for computational simplicity, one representative τ 0 is chosen. A simple method is to choose the same τ for all source-detector pairs [52]. We use a τ 0 value in which g 1,c (r s ,r d ,τ) decays to e −1 for each source-detector pair, since it yields the optimal data sets [45]. The steps to solve the inverse problem with a non-linear solver are illustrated in Fig. 2. First, a finite element mesh is utilized to define the geometry. Then, an initial guess of the flow distribution is made and the optimal data set is selected. The iteration to update flow estimate starts with the calculation of g 1,c (r s ,r d ,τ) based on the flow estimate and is followed by obtaining the difference between g 1,c and g 1,m , φ i = ln(g 1,m (r si ,r di ,τ) / g 1,c (r s ,r d ,τ)). After that, the stopping criteria are evaluated. The stopping criteria are set such that when the decrease in χ 2 is less than 2% for two successive reconstructions, the iteration stops. However, the final result is not selected as the result from the last iteration. More detailed discussions on selection of optimal results using L-curve analysis are presented at the end of this section. If the stopping criteria are not met, a Jacobian matrix, W, is built, which satisfies is the difference between flow estimate and true flow distribution. Note that W is built in the same way as been described in our previous where T indicates a matrix transpose and η is a regularization parameter which is decreased with each iteration. The formulation includes Tikhonov regularization to stabilize matrix inversion [53], which is especially effective for solving the ill-posed equations (i.e., the number of unknowns is much larger than that of linear independent equations). As the final step in each iteration, Δ(αD b (r)) (n) is added to the flow estimate, yielding αD b (r) (n) . The iterative update of the flow is carried out until the stopping criteria are met. An inverse problem solver generally needs to minimize χ 2 . However, with the existence of noise in the measured data and the ill-posedness of the inverse problem, iterative minimization of χ 2 may amplify noise into image artifacts such as a flow estimate with unreasonably high values (e.g., more than 50 times the background flow) at some positions. Thus, the final iteration of the iterative reconstructions may not necessarily provide physiologically reasonable flow images although it gives the smallest χ 2 . To identify the occurrence of this problem and to choose the optimal reconstruction results, an L-curve analysis is performed and served as the final criteria to select the optimal reconstruction result [54].
With each iteration in the reconstruction, a total flow perturbation from the initial guess can be calculated as δ(αD b ) (n) = αD b (r) (n) -αD b (r) (0) . The norm ||δ(αD b ) (n) || provides a measure of the fluctuation in the reconstructed flow estimate. On the other hand, the χ 2(n) calculated at each iteration, as defined in Eq. (4), measures the difference between g 1,c and g 1,m and quantifies the quality of the flow estimate αD b (r) (n) . For the L-curve analysis, the cost function χ 2(n) at each iteration is plotted against the norm of the total flow perturbation ||δ(αD b ) (n) ||, which gives an "L" shape curve. To account for the different orders of magnitude between the two axes, both parameters are normalized to their respective maximum values. Then, the curvature of the L-curve at each data point is calculated. Finally, the iteration which generates the data point with the maximum curvature in the "L" curve is selected as the optimal iteration and the corresponding αD b (r) is taken as the final reconstructed flow distribution. An example of the typical L-curve from reconstruction on mouse data is shown in Fig. 3. Fig. 3. A L-curve obtained from the iterative image reconstructions of mouse data. L-curve analysis is employed to choose the optimal results from image reconstructions, which is the iteration with the maximum curvature on the L-curve as shown in red dot.

Numerical simulation to compare a linear and a non-linear reconstruction method
To compare the non-linear inverse problem solver with a linear solver that we reported previously [47], reconstructions based on simulated data were carried out. g 1 of multiple source-detector pairs over a region of heterogeneous flow distribution was calculated using a FEM based solver of Eq. (2) modified from Near Infrared Fluorescence and Spectral Tomography (NIRFAST), a software package for diffuse optical tomography [55]. The simulation geometry consists of a cylindrical region with higher flow and a surrounding homogeneous background, as shown in Fig. 4(a). It is discretized using a finite element mesh (x: −8.0 mm to 8.0 mm, y: −13.5 mm to 13.5 mm, z: 0 mm to 8.0 mm, total number of nodes: 201791, total number of elements: 1031591). The higher flow region has a diameter of 1.5 mm and a length of 15 mm, and it is submerged 2 mm deep into the background. The cylindrical region has αD b = 2.0x10 −6 mm 2 /s, and the background region has αD b = 1.0x10 −6 mm 2 /s. Figure 4(b) shows the pattern of the source and detectors utilized to scan over the region for data collection, which has source-detector separations of 3.0, 4.5, 6.0 and 7.5 mm. The pattern was scanned in both x and y directions, with step sizes Δx = 0.35 mm and Δy = 0.70 mm. The scan coverage is shown in Fig. 4(c), which only includes source and detector positions inside the rectangular region of −3.5 mm < x < 3.5 mm and −10 mm < y < 10 mm in order to speed up the calculation. The original pattern consists of 1764 source-detector pairs and the final pattern includes 1272 pairs.
To compare the reconstruction results, a coarse mesh (x: −8.0 mm to 8.0 mm, y: −13.5 mm to 13.5 mm, z: 0 mm to 8.0 mm, total number of nodes: 94826, total number of elements: 462047) was used for both the linear reconstruction method and the non-linear reconstruction method. No noise was added to the simulated data.

In vivo mouse measurement protocol
For the in vivo bone graft measurements, three groups of mice were used for DCT imaging: the autograft group, the allograft group and the T.E. allograft group. The animal study was conducted according to a protocol approved by the University Committee on Animal Resources (UCAR) at the University of Rochester. The details of the graft surgical procedures can be found in the literature [56, 57]. Briefly, a 4 mm segment of the left femur of a BABL/c mouse is removed and replaced with a graft. For the autograft, the removed segment was transplanted back. For the allograft group, decellularized bone segments from C57BL/6 mice (different strain than BALB/c) were utilized for the grafts. For the T.E. allograft group, decellularized bone segments from C57BL/6 mice were seeded with MSCs using hydrogel [57], and then transplanted into the femoral defect. For each group, 5 mice were initially planned for the study, since the purpose of the study was to identify the spatiotemporal characteristics of the blood flow index of different groups and provide information needed to calculate sample size for future studies. However, one mouse in the allograft group and one mouse in the autograft group died due to the adverse effects associated with ketamine, anesthetic for the surgery, before the start of any post-surgery DCT measurements, and thus were excluded from the study.
As illustrated in Fig. 5(a), DCT scans were performed before the graft surgery (denoted as week 0), and then every week after surgery for 9 weeks. For each scan, the mouse was placed on a warming pad and anesthetized with continuous isoflurane administration (1 -2.5%). Isoflurane was chosen for minimal interference with microcirculation, fast anesthetic effect and reliability [58-60]. Fur was removed from the scanning site using electric clippers. Then, the source-detector pattern shown in Fig. 4(b) was scanned over the tissue surface with a step size of 0.35 mm in the direction perpendicular to the length of the femur and 0.70 mm in the direction parallel to the length of the femur. The scan region covered the entire femur with the graft, with the femur located at the center, similar to the simulation geometry shown in Fig.  4(c). The four corner positions of the rectangular scan area covered by the source were marked using a skin-compatible surgical marker during the first measurements (i.e., presurgical measurements) and maintained throughout the longitudinal measurements. Alignment of the scanning pattern with the co-registration markers was performed before each DCT scan thereafter. At week 11 after surgery, an additional DCT scan was performed and then followed by micro-CT imaging of the mouse leg at the Center for Musculoskeletal Research at the University of Rochester. To co-register DCT scans with micro-CT, barium sulfate markers were applied on the mouse skin on top of the DCT scan markers and served as reference points in the micro-CT images. After the micro-CT scan, standard Digital Imaging and Communications in Medicine (DICOM) format medical images of left hindlimbs were obtained.

Generation of FEM mesh based on anatomical information from micro-CT
The micro-CT images were processed with NIRFAST Slicer [61], a free software package which can view medical images, perform de-noising and image-segmentation, and export segmented medical images into Matlab (Mathworks, Natick, MA) for 3D finite element mesh generation. After the images were imported into NIRFAST Slicer, de-noise algorithm with median filter was performed and then the images were segmented using a built-in K-means and Markov random field algorithm [61], and divided into the region of bone, soft tissue and air. The segmentation was fine-tuned manually using the paintbrush tool. The reference points for the DCT scan region were marked on the segmented images. Finally, segmentation results were sent to Matlab, by which 3D finite meshes of the legs were generated. The generated mesh had the region of the transplanted graft, the remaining host femur and the surrounding soft tissue labeled as different regions and typically consisted of more than 100,000 nodes in total, as is shown in Fig. 5(b).
After the 3D mesh was generated in Matlab, the source-detector pattern from DCT scan was registered on the mesh by overlaying the co-registration markers with the corresponding source positions and linearly interpolating the other sources and detectors on the tissue surface. The tissue surface curvature was not accounted for in current study as it was usually small. In addition, the sources and detectors within 1 mm of the mesh boundaries, where the tissue curvature was big, were excluded. The same absorption coefficient μ a = 0.01 mm

Statistical analysis
To determine the existence of statistical significance, Shapiro-Wilk test was utilized to test the normality of the data distribution for each group. One-way ANOVA or non-parametric Kruskal-Wallis test was utilized based on the result of the normality test. When the null hypothesis (i.e., three groups have the same mean) was rejected (p < 0.05), Tukey's honestly significant difference procedure was used for multiple comparisons to identify the pair of groups which demonstrated statistical significance.

Comparison of a linear and a non-linear reconstruction method
Reconstruction results from a linear DCT reconstruction method and a non-linear reconstruction method are shown in Fig. 6. To compare the results, we define the region that has the final flow change larger than 1/e of the maximum change as the recovered anomaly region. The expected center of the anomaly region is (0.0, 0.0, 2.0) mm. For the linear solver, the reconstructed anomaly region centers at (−0.2, −0.5, 1.4) mm; for the non-linear solver, the recovered anomaly region centers at (0.2, −0.9, 1.8) mm. Compared with the linear solver, the non-linear solver gives a much better depth localization of the reconstructed anomaly. However, results from the non-linear reconstruction still need improvement. For example, the average αD b value within the anomaly from the non-linear reconstruction is lower than the expected value and the high αD b region disperses beyond the expected z range (1.25 -2.75 mm).

Relative blood flow changes in a mouse with an autograft
The longitudinal changes in 3D rBF of a mouse with an autograft were shown in Fig. 7 as an example of DCT reconstructed flow changes. In Fig. 7, the rows show the rBF distribution in weeks 0, 1, 3, 5 and 7. The columns show the rBF at the tissue depths equal to 0.5, 1.5, 2.5, 3.5, 4.5 and 5.5 mm. The profile of the bones, including the femur, the tibia and the coxal bone are outlined using black lines and the graft region is shown with red lines. The same colorbar was used for all plots. From the results, it can be seen that the blood flow in an autograft mouse peaked at week 1 and gradually decreased as healing progressed. At the final stage of healing, blood flow returned to the pre-transplantation baseline level. Moreover, blood flow increase was mainly localized in the graft region for the autograft mouse, as shown in the tissue depth of 3.5 mm. At shallow depths (as shown in slices for z ≤ 1.5 mm), a few regions of high blood flow were observed in the soft tissue region. On the contrary, at deeper region with the graft, the changes in soft tissue blood flow were much smaller compared to those in the graft region (as shown in slices for z ≥ 3.5 mm).

Comparison of blood flow changes between different groups
To compare the difference in blood flow changes among the three different groups, images from one representative mouse in each group at the tissue depth passing through the center of the femur were shown in Fig. 8. The three rows show the results from the autograft group, the allograft group and the T.E. allograft group, respectively, and the columns show different time points at weeks 0, 1, 3, 5 and 7, respectively. The outlines of the bones and the grafts are shown in black and red lines, respectively. The same colorbar was used for all plots.
Compared to the blood flow distribution in the autograft mouse (Fig. 8 top) that showed a well-localized blood flow elevation within the bone primarily in the graft region at week 1 and 3, blood flow in the allograft mouse (Fig. 8 middle) exhibited multiple regions in soft tissues with higher blood flow than that in the graft. At week 5 and 7, blood flow in these multiple regions decreased and more localized blood flow elevation was found within the allograft.
The blood flow in the T.E. allograft mouse (Fig. 8 bottom) exhibited spatial and temporal features similar to both autografts and allografts. At week 1, blood flow elevations were found within the graft as well as in multiple regions in soft tissue. The level of blood flow elevations was similar between the graft and the soft tissue at week 1, but the elevation within the graft became more prominent compared to soft tissues as healing progressed.

Average rBF changes in graft region for different groups
The longitudinal changes of the average blood flow in the graft region for an individual mouse were shown in top row of Fig. 9. As mentioned previously in Section 2.4, there were 4 mice in the autograft group, 4 in the allograft group and 5 in the T.E. allograft group.
In the autograft group, graft blood flow increased significantly from the baseline level 1 week after the graft transplantation and decreased more than half of the initial elevation at the 2-week time point. After that, blood flow decreased gradually towards the baseline level. This pattern was consistent in every autograft mouse. For the allograft group, more variations between the mice within the group were observed. The time to peak ranged from week 1 to week 5. The average graft rBF returned to the baseline at a slower pace than that of the autograft group. For the T.E allograft group, the blood flow peaked at week 1, similar to that of the autograft. However, the trend in blood flow decrease was not as consistent as that of the autograft group. Overall, the changes in the T.E. allograft group were more similar to the autograft group compared to the allograft group in that 4 out of 5 mice demonstrated a similar temporal pattern to the autograft.
The bottom row of Fig. 9 shows the group averaged rBF within the graft. Even though different groups showed varying degrees of averaged rBF at week 1, it was not statistically significant partly due to the limited sample size. On the other hand, observations made for individual mouse strongly suggest that whether the rBF peaks at week 1 can be a potential criterion to distinguish the degree of healing. To quantify this, we defined Percentage of Peak (POP) at week 1 for each mouse, as the following: where max(rBF(t)) is the maximum rBF from the longitudinal curve spanning the monitoring period as shown in the top row of Fig. 9. The group averaged POP at week 1 are shown in Fig. 10, along with the standard error of the mean as error bars. The statistical test on POP revealed that the autograft group is significantly different from the allograft group, but not from the T.E. allograft group.

Discussion
In this study, we used an improved version of the DCT system from our previous publication [47], to monitor the longitudinal blood flow changes in mice with femoral grafts. Using the non-contact scanning system and a non-linear reconstruction method, we showed spatial and temporal differences in blood flow among three groups of mice with different graft types over 11 weeks of healing. The study is unique in that: (1) the non-invasive nature of the DCT instrument enabled frequent longitudinal blood flow monitoring of the same mouse throughout the experiment; (2) DCT based on spatially dense data and non-linear reconstruction algorithm provides comprehensive spatial and temporal blood flow changes throughout the healing process, which was not available in previous studies due to limitations of other imaging modalities; (3) the study revealed different spatial and temporal blood flow changes in autografts, allografts and T.E. allografts with MSCs, strongly suggesting blood flow changes in early stages of healing (< 5 weeks) may be utilized to predict bone healing outcome of different graft types.

Longitudinal blood flow changes in grafts
The longitudinal rBF changes show a general trend of blood flow elevation after the graft transplantation and subsequent return to the pre-surgery baseline level. However, the autografts, the allografts and the T.E. allografts show different temporal patterns. The autograft group consistently displayed a peak in graft blood flow at 1 week after the transplantation. From the tomographical results, the flow changes in autografts were also well localized in the graft region. The allograft group, on the other hand, exhibited more variations in the change of blood flow within the group, with blood flow peaking from 1 week to 5 weeks after the transplantation and the elevation returned to the baseline at a much slower rate. In the 3D results, more elevations in the tissue surrounding the graft were observed compared to the autograft group. The T.E. allograft group showed intermediate changes between the autograft group and the allograft group, with blood flow peaking at week 1, similar to the autograft group: however, the rate of blood flow decrease was slower. The blood flow changes were better localized at the graft region compared to the allograft group. The trend in blood flow changes in autografts, which showed a sharp blood flow increase followed by a return to the baseline level, agrees with previous studies in bone graft/fracture blood flow monitoring, including studies using radio-active microsphere [63], laser Doppler flowmetry [18], and DCS [56]. The increase in blood flow at week 1 is likely due to several factors including inflammatory response to surgical insult, hyperemia, and angiogenesis associated with tissue regeneration [64]. It is worthwhile to note that, since our measurements were conducted on a weekly basis, the immediate drop of blood flow (within 3 days after the surgery) in the graft following the surgery reported by other techniques was not captured in the current study. In addition, the observation that the allograft group exhibited a longer period of blood flow elevation is consistent with previous results observed by DCS [56]. The elongated blood flow elevation may be an indicator for delayed healing in allografts, which has been observed in other studies evaluating the bone callus formation and vasculature formation [13,57].
With the rich spatial and temporal information provided by the non-invasive DCT, differences in blood flow changes in autografts, allografts, and T.E. allografts with MSCs were clearly captured, as shown in Fig. 9. Autografts exhibited consistent blood flow elevation at week 1 which subsequently decreases, consistent with good healing. Although not identical to autografts, the T.E. allograft group showed a similar blood flow pattern to autografts compared to the untreated allograft group, demonstrating improvements in bone graft healing using this tissue engineering approach. Hoffman et al. demonstrated that a mixture of unaltered MSCs and osteoprogenitor cells can further improve the bone callus formation, bridging endochondral bone formation and bone bio-mechanical stability [13]. In the future, DCT may play a valuable role in assessing the healing potential of T.E. allografts with different compositions of cell populations, soluble cues, and/or matrices to guide the optimization of tissue engineering approaches.

Development of methods
In the methods and the results sections, the improvement on instrumentation and image reconstruction method are presented, which includes: (1) increased light throughput in the non-contact scanning system, enabling the collection of spatially dense data at multiple source-detector separations (up to 7.5 mm); (2) a non-linear reconstruction method that gives better depth localization, as opposed to the previous linear reconstruction program; (3) use of micro-CT derived 3D mesh and assignment of optical properties based on anatomical information, compared to semi-infinite mesh assumption used previously [47]. With the current methods, longitudinal blood flow changes during graft healing were successfully monitored non-invasively and the spatial and temporal differences among three graft types were revealed. The versatility of the current system also enables direct application in other studies, such as monitoring tumor blood flow change after chemotherapy in rats [65]. Despite the improvements, several limitations remain to be addressed in future studies.
First, the non-linear reconstruction method has limitations. The simulation results show image artifacts in the shallow depth (<1.5 mm) and error in the reconstruction position and contrast. To suppress the image artifacts, spatially variant regularization to account for the sensitivity difference at different tissue depths might be used [45]. The reconstruction may be further improved by better utilization of the anatomical information. To this end, reconstruction with soft-priors will be implemented in the future as it has been shown to be more robust to noise and give better reconstruction results [66]. In addition, in the average rBF results (Fig. 9), mouse 7 exhibited much smaller changes compared to other mice. Inspection of the anatomical structure revealed that the femur of mouse 7 (centered at z = 5.5 mm) was much deeper than the others (typically centered at z = 3.5 mm). Considering the maximum source-detector separation used is 7.5 mm, the changes in the graft may not be detected as well as other mice with shallower bones. This may imply the need for a more rigorous study on the depth sensitivity of the current system and the necessity of criteria for correction or exclusion based on the bone depth.
Second, the influence of the tissue surface curvature on the non-contact data collection was not accounted for. The non-contact scanning system projects a planar pattern for data collection. However, in the image co-registration, the source and detectors are interpolated to the tissue surface and their heights followed the profile of the tissue, which is not planar. To partially account for this inconsistency, the sources and detectors within 1 mm of the mesh boundaries, where the height changes are big, were excluded from the reconstruction. Despite the exclusion of those data, the influences of tissue curvature on image reconstruction were not systematically investigated in the current study. A potential solution is to account for the propagation of the auto-correlation function in free space, similar to a method proposed for diffuse optical tomography [67].
Third, the image co-registration scheme is not ideal. To record the DCT scan region and co-register it with the micro-CT images, DCT scan markers and barium sulfate markers were used. The size of the markers, especially the Barium Sulfate markers, which are typically 1 mm in diameter, is comparable to the scanning step size. Although only the center of each marker was used, possible mismatch is expected. To study the effect of the mismatch in image co-registration, a reconstruction on the autograft mouse data as shown in Fig. 7 with a 0.5 mm displacement of the source-detector pattern was carried out. The results showed the displacement in source-detector pattern mainly introduced a spatial shift of the high blood flow region for about 0.5 mm in the reconstructed results. The general distribution and the temporal trends did not change. In the future, a more accurate image co-registration scheme may be needed and a more rigorous study needs to be designed to study the stability of the reconstruction with image co-registration mismatch.
Finally, the assignment of optical properties from literature is not ideal. Optical properties, especially reduced scattering coefficients μ ś can greatly influence the accuracy of reconstructed αD b values [47]. To account for the difference in scattering in soft tissue and bone, μ ś from literature was assigned in the current study. In the future, better optical property assignment may be achieved by combining the current DCT system with a diffuse optical tomography system that can provide 3D absorption coefficients and reduced scattering coefficients of each mouse.
Additionally, there is a potential concern related to the use of the Siegert relation in our analysis. It has been reported in the literature that the signals from adult human bone may not satisfy all the assumptions of the Siegert relation under certain conditions [68]. In our case, the necessary conditions to satisfy the assumptions are provided by the presence of a thick muscle layer over the bone and natural ensemble averaging through breathing-induced motion [22,56]. In the future, the limits of the approximations in small animal models could be explicitly checked.
In spite of the limitations, the current DCT study provided previously unavailable information on the 3D longitudinal blood flow changes in three different graft types. Recently, there has been an increased research interest in the development of DCT scanning system and application of DCT in different studies. He et al. used a motorized stage to move the source and detectors over the scan region and used this system to capture the blood flow contrast in patients with breast tumors [46]. Johansson et al. constructed a scanning system that used a motorized stage to move the subject rather than the source-detector pattern and combined diffuse correlation measurements with broadband diffuse optical measurements to study the hemodynamics of mouse tumor [69]. However, in general, DCT is still a technology that has not been widely applied. To this end, we plan to compare our DCT measurements with more established methods, such as fluorescent microspheres or dynamic contrast enhanced MRI in the future.

Summary
In summary, we have utilized a DCT system to study the longitudinal blood flow changes in three groups of mice with different graft types; i.e., autografts, allografts and T.E. allografts. Our system enabled the collection of spatially dense data with multiple source-detector separations for 3D imaging reconstruction, utilized a non-linear reconstruction approach to improve the accuracy in the reconstruction and leveraged the anatomical structure of each mouse obtained with micro-CT to define the geometry and optical properties. The DCT results show that after the graft surgery, blood flow increased in grafts for all groups and then returned to the pre-transplantation baseline level. The autograft group exhibited consistent trends in blood flow changes in the graft region, with the blood flow peaking 1 week after surgery. The allografts showed more temporal variations in blood flow within the group. The T.E. allografts exhibited blood flow changes that were more similar to the autografts. The results monitored by our DCT system strongly suggest that 3D longitudinal blood flow changes may serve as surrogate markers for assessing the success of bone healing at relatively early time points.