Biological tissue both absorbs and scatters light. Diseased tissue can be distinguished from normal tissue on the basis of differences in optical properties arising from physiological variations. Optical imaging with probes formed by biologically relevant molecules targeted at an intrinsic molecule or receptor overexpressed in a disease conjugated with a fluorescent marker to visualize the receptor can provide a cost-effective molecular imaging mechanism. The ability to image and visualize molecular scale events with targeted probes makes optical imaging an attractive modality for small animal imaging. Non-ionizing light facilitates the use of optical probes for longitudinal in vivo studies and provides an effective platform for drug discovery and development.1, 2, 3, 4
The detectability of fluorescent contrast, image reconstruction, and ability to visualize and distinguish between normal and diseased tissue is influenced by the nature of light transport and distribution. Molecular imaging necessitates careful design of optical instrumentation and application of appropriate fluorescent probes in adequate concentration to enable the detection of pertinent information from the scattered signal.
In this paper, we present a generic method based on optical transport to create a digital mouse phantom and illustrate its usefulness in predicting the performance measures for an imaging system. A phantom could be used to design imaging systems, assess imageability of fluorescent probes, and aid in model-based iterative optical image reconstruction. Predictive design could also help in preventing the sacrifice of animals and in bridging the gap between animal models and clinical studies.
A digital phantom is a numerical representation of a structure that can be represented geometrically. Advances in computing technologies over the past two decades have replaced the early analog phantoms and have led to a spurt of digital phantom–based simulation of imaging systems.5, 6, 7 Phantoms have been used extensively in designing model-based image reconstruction algorithms for optical imaging. Proofs of concept experiments have involved the use of both canonical geometries and finite element (FE) meshes to model anatomical structures. For example, phantoms ranging from finite rectangular slabs,8, 9 infinite slabs,10 circular or cylindrical geometries,11 and three-dimensional (3-D) FE mesh for a conical shape12 have been used to model highly scattering properties of breast tissue.
Magnetic resonance imaging (MRI) images have been used as digital priors for over a decade.13, 14 Use of prior anatomical information in the form of MRI images to improve optical reconstruction was first suggested by Barbour 15 Chang used a linear perturbative model for time-independent optical sources for MRI-assisted optical reconstruction for mammography.16 MRI priors have been recently proposed for the identification of probability density functions that are used in a formulation of optical reconstruction within a Bayesian framework.17
Optical reconstruction for brain imaging for the rat cranium with MRI priors was reported by Pogue 18 The approach involved the use of a fine FE mesh of the segmented MRI image with prescribed optical properties to perform optical measurements and a coarse FE mesh for the reconstruction. MRI brain image slices segmented into regions depicting skin, bone, and gray and white matter have been used to create a layered head model with FEM meshes representing anatomical priors for use in optical reconstruction.19 3-D digital maps of an MRI rat cranium have been used to perform hemodynamic studies of the rat brain.20
In a different application involving the visualization of subsurface lesions with near infrared (NIR) cameras, FE digital phantoms from MRI images have also been used in novel volume rendering techniques that mimic NIR light transport in a human arm.21
Adaptive FE meshes that incorporate varying mesh resolutions based on the shape and size of anatomical structures and targets have been pioneered in model-based iterative optical image reconstruction schemes by a number of researchers.22, 23, 24, 25
In a recent publication, Dogdas 26 have reviewed the existing small animal imaging phantoms and described the process of creation and design of Digimouse, a digital mouse phantom generated from co-registered x-ray CT and cryosection data, and applied it to simulation of bioluminescence.
A diffusion approximation (DA) to the radiative transfer equation (RTE) has been extensively proposed and used to model photon transport in biological media. Analytical expressions of DA for infinite cylinder and slab geometries in the time domain have been formulated and applied to optical image reconstruction. FE equations of DA in the time and frequency domain in two dimensions and three dimensions have been presented and utilized in absorption, fluorescence reporter concentration, and lifetime reconstruction as well as in small animal scanning systems and in multimodality imaging with ultrasound and MRI for clinical applications in brain and breast imaging. 15, 16, 17, 18, 19, 27, 28, 29, 30, 31, 32, 33, 34
DA is known to be applicable when scattering is much greater than absorption and at distances greater than a few scattering lengths from the source. In situations pertaining to low scattering regions and subsurface imaging, radiative transfer equation and a stochastic Monte Carlo model for light propagation in tissue are the more exact formulations. 35, 36 However, simplicity, ease of use in the case of complex geometries, and potentially faster computation are the key reasons why DA continues to be the widely preferred biomedical photon transport model.
We utilize our validated numerical implementation of the 3-D diffusion equation to design the digital mouse phantom.37, 38 We first create the boundary of the mouse phantom from MRI data using a series of image processing steps. We next discretize the geometry with an FE descriptor using FEMLAB. The mesh resolution of the FE model is representative of the optical properties of biological tissue and optimized to perform on a computing system with RAM. We then demonstrate through simulations that a section of the mouse near the light source is adequate to predict performance of an imaging system and that the variation in intensity of light detected on the boundary is well within typical noise levels for up to 20% variation in optical properties and number of nodes used to model the boundary of the phantom. Last we illustrate the significance of modeling the undulating boundary of a mouse by predicting a tenfold change in the detectability of a -deep fluorescent target dependent on the local curvature of the boundary. The range of detectability has implications for the desirable specific binding of targeted fluorescent contrast agents based on mouse anatomy. This work is aimed at standardization of digital phantoms used in optical imaging system development. To the best of our knowledge, this is the first study that utilizes a comparative analysis of the performance of 3-D FE diffusion and 3-D analytical diffusion models on a cylindrical geometry to assess the sensitivity and use of fluence computations on an FE model with changes in the optical properties and boundary of a small animal model.
Our methodology has been inspired by the extensive use of MRI data as anatomical priors and the diffusion equation as a model for photon transport in biological tissue. The broad steps consisting of generating mouse geometry from MRI data, using predictive fluence computations to arrive at optimal size and mesh resolution of the phantom, and assessing its robustness with changes in optical properties and boundary are shown in Fig. 1 . Figure 1 also depicts the principle sources of variability likely to affect the final design of the phantom in each step.
Creation of Mouse Phantom
We have generated the mouse phantom from MRI data using the following procedures.
We extract the mouse boundary from MRI data through the application of image processing algorithms, as shown in Fig. 2 . The MRI data is first smoothed using a one-dimensional (1-D) Gaussian filter. A threshold is applied next on the filtered image volume. The crispness of the extracted boundary is particularly sensitive to the value of the threshold. Lower threshold values result in a boundary cluttered with undesirable clusters of points, whereas the higher threshold value results in the extraction of internal organs of the mouse along with the final boundary. In most real-world problems, foreground and background regions do not follow perfect bimodal intensity distributions, resulting in some overlap of intensities. We choose a threshold that eliminates as much of the background as possible without affecting the foreground object. Noise tends to persist despite thresholding and is subsequently removed through a process of erosion, image subtraction, and connected component analysis.
The original MRI data set of the mouse has 26 slices, each with dimensions . Resolution within each 2-D slice is . The adjacent slices are separated by . We do not consider some of the frames toward the two edges of the mouse for geometry creation because the entire boundary is not retrieved in these slices. The geometrical representation of the mouse is created from the extracted boundary through a process of selecting frames, contours, and edge segments and finally lofting the frames to generate a 3-D structure. For every two-dimensional (2-D) frame, we first generate a contour followed by a 2-D curve and a 2-D solid mesh. Last, we loft each 2-D solid mesh with appropriate spacing to create the geometry, as shown in Fig. 3 .
We use FEMLAB39 to create the geometry and FE mesh of the mouse phantom. The transformation of mouse geometry to mouse mesh is depicted in Fig. 3. Mesh resolution is dictated by the needs of the diffusion model to accurately simulate light propagation in biological tissue and constrained by the available hardware and software used for FE analysis. Our earlier studies to map the validity of using an FE 3-D diffusion model as an engineering tool for the predictive assessment of imageability through a systematic comparison of depth-resolved fluence calculated by diffusion equation and the gold standard Monte Carlo model had shown that although, in theory, the diffusion model applies when , in practice, its performance can vary sharply for the same ratio.37, 38
We summarized the rules of thumb for applicability of the FE diffusion equation for excitation fluence as:
1. The average mesh resolution should be higher than the mean free path close to the source in the region of interest (the average value of node-to-node distance ).37
3. The model is not applicable37 at distances .
4. The model shows the maximum departure37 from the Monte Carlo at a distance of . This difference can be estimated so that appropriate tolerances may be used in design. The mean free path and penetration depth are given by:is the absorption coefficient, is the scattering coefficient, is the anisotropy, and the reduced scattering coefficient .
Phantom Size and Mesh Resolution
Range of optical properties in biological tissue.
|μs (mm−1)||g||μs′ (mm−1)||μa (mm−1)||nr|
We first construct an equivalent cylindrical phantom with the same length and similar mesh resolution as that of the desired mouse phantom. In accordance with the rules of thumb, a single cylindrical phantom capable of simulating light propagation across the range of optical properties would need to have the size specified by the highest penetration depth and the mesh resolution specified by the lowest mean free path. This phantom should have at least a diameter of and a length of , and its FE model should have an average node-to-node distance of . Our current hardware ( HP UNIX workstation) can provide a maximum mesh resolution of for a cylinder of this size. Based on the capability of the hardware, we have designed a cylindrical phantom of diam and height represented by 199,414 nodes, 1,146,319 elements, and 30,690 boundary nodes.
The digital mouse phantom shown in Fig. 4 is constructed by simply shortening the mouse geometry by an equal distance from the axial extremities for the mouse phantom to have the same length and mesh resolution as that of the cylindrical phantom.
Based on our earlier studies comparing light fluence computed by Monte Carlo simulation and the diffusion equation, we expect the attributes of the “design” phantom to be appropriate for prediction of photon transport for , , and . In order to verify and assess the applicability of the phantom for light simulation at varying axial and angular distances from a point source located on the boundary, we measure simulated fluence on the boundary of five slices axially located at and shifted by , on either side of a point source positioned approximately near the midaxis.
We use the comparative performance40 of the FE diffusion equation and our MATLAB implementation of a 3-D analytical diffusion model27 on a cylindrical geometry to verify the size and mesh resolution of the mouse phantom for the “design” optical parameters. We evaluate excitation fluence from the FE model (Fe) and compare it with the analytical solution (Fa) at all boundary points in each slice for the optical properties corresponding to the phantom by computing where is the axial coordinate and is the source-detector angle in the -plane. We then perform a sensitivity analysis of the “design” phantom by increasing the scattering and absorption coefficient by 5, 10, and 20% to study the variation in . The results of the sensitivity analysis are presented in the next section.
As shown in Fig. 1, there are several image processing steps, each of which is a potential source of variation in the generation of a mouse phantom. For the current study, we focus on the number of edge segments that are used to describe the contours selected as the boundary for each frame. Local geometry and mesh resolution close to the source and detectors are known to have a strong influence on simulated fluence.11, 38 Sources and detectors are located on the boundary. Therefore, it is pertinent to assess the sensitivity of the simulated results to variations in the boundary modeled through a change in the number of edge segments. We consider the number of edge segments in the basic mouse design and generate mouse meshes with edge segments such that the number of boundary nodes change by approximately , , and from the number of boundary nodes in the designed case. The negative variation is generated by down-sampling the edge segments from the base design, and positive variation by up-sampling the edge variation from the base. The base design has 60 edge segments. The change in the number of segments leads to a change in the number of nodes in the FE model. The number of edge segments, nodes, and segment lengths in each design and their percent difference from the base design is shown in Table 2 .
Attributes of design variants of the mouse model.
|Design||No. of Edges||Nodes||Boundary Nodes||Segment Length|
|4||80||133300 (1.1%)||30814 (6%)||0.9752|
|5||95||134596 (2.1%)||32360 (11%)||0.9720|
|6||100||138389 (5%)||34996 (20.4%)||0.9532|
We compare the excitation light fluence for the six variants in boundary models with the base design in five slices axially located at and shifted by , on either side of the point source positioned approximately near the midaxis, similar to that in the cylindrical phantom. As shown in Fig. 4, the source is located in the central slice 3. Each slice has a thickness of . The results of the comparison are presented in the next section.
Verification of Design Parameters
We compute the variation in excitation fluence40 with source-detector angle from our 3-D FE implementation of the diffusion model (Fe) and the analytical solution of the 3-D diffusion equation (Fa) on the cylindrical geometry in all five slices over a range in the source-detector angles for the “design” optical properties ( , , and ), ensure that at distances that are a few away from the source,11 and observe that as the distance from the source , .
Sensitivity to Variation in Optical Properties
For our analysis, we eliminate all nodes in the slice on the boundary that are within 2 of the source and compute mean ratios for sixteen representative cases over the five slices with mean free paths corresponding to 0%, 5%, 10%, and 20% increases in and from their “design” values of and , respectively. The mean ratio is found to vary linearly with over the design space. The mean ratios increase with increase in scattering coefficient and decrease with increase in absorption within each slice.
The differences in the FE and the analytical diffusion model could arise from both meshing and the differences in the manner in which the boundary condition is incorporated in them.3, 27 We would expect to represent the average pattern of comparison within an error margin. The error is defined as:and are the mean ratio of the “design” case and the test cases, respectively, calculated over all the slices. The error for each test case used in the sensitivity analysis is shown in Table 3 .
The error rises to a maximum of 5.7% for up to a 20% increase in the optical absorption and reduced scattering coefficients for source-detector distances varying from to about . However, as increases to ( increase), our computations (not presented here) show that at a source-detector angle of , FE fluence reduces seven orders of magnitude faster than analytical fluence at , four orders of magnitude faster at , and two magnitude faster at . It could be expected that at large distances from the source, the FE fluence would converge to analytical fluence predictions and this trend would prevail as long as the photon transport is diffusive, i.e., . In measurements of fluorescence during frequency domain photon migration, Thomson and Sevick-Muraca have reported a mean accuracy of 5.4% in modulation depth and in phase up to distance from a point source. Variability in measurements increases to 17% in modulation depth and in phase for longer distances from the source.41 Our simulated variations in excitation fluence are of the order or less than typical measurement noise and permit the use of the cylinder size and mesh resolution to design the mouse phantom for optical properties within 20% of the design value.
Sensitivity to Variation in the Boundary
Since there is no point-to-point correspondence of boundary nodes among the six variant mouse phantom designs that we described in the preceding section, we subdivide each slice into four quadrants for the purposes of comparison. We use the mean of log of fluence in each quadrant in each slice as a representative estimate in each design and compare those estimates with the same measure in the basic design. Thus, the basis of our comparison is a ratio of the mean of log. Mean of log as an estimate makes sense since the fluence varies exponentially with source-detector angle.
In order to understand the nature of the variation, if any, between each design variant, we study the similarities and/or differences across the slices of each design. We ensure42 that each of our data sets follow a normal distribution while rejecting outliers that lie within 2 of the source. We use analysis of variance (ANOVA) tests42 to assess whether there are significant differences of ratios of mean of logs over all slices and all variant designs, over each slice across all designs, and over each design across all slices. The null hypothesis for the two-way ANOVA test is that the ratios of mean of logs over all mouse designs in all regions do not differ significantly. At a 95% confidence interval, we find no statistically significant difference between our six design variants. Change in light intensity as we move away from the source in each design is a more dominant source of variation compared to the change in the FE model of the mouse . This is also evinced by our results (not shown here) of a one-way unstacked ANOVA test for comparing the similarity of distribution of ratios within designs and within slices.
One of the important design goals for exogenous fluorescent contrast agents is to achieve an adequate specific binding to the disease target to enable detection through an optical imaging system.43 Target-to-background ratio (TBR) is a measure of specific binding. TBR is a function of the sensitivity of the imaging system, optical properties in the interior of the tissue at the wavelengths of excitation and emission of the fluorescence, time after injection of contrast agent, and depth and size of the anomalous tissue and could be simulated using tissue mimicking phantoms.
To illustrate this, we embed a cylindrical fluorescent target of height and diam, beneath the boundary of both the cylindrical and the mouse phantoms and illuminate it with an excitation point source placed on the boundary above the target. The absorption coefficient of the target is tenfold higher than the background and is assumed to be the same at the excitation and the emission wavelengths. We use our 3-D FE diffusion model to calculate both excitation and emission fluence in the reflectance geometry on the boundary close to the source for the two different target positions. The local geometry near the source for the two positions is shown in Fig. 5a . The heights are averaged over the nodes that are within an angular distance of on the boundary over every axial distance for on either side of the source (at for purposes of comparison). A similar averaging is used to plot the excitation and emission fluence profiles, shown in Figs. 5b and 5c.
Percent error for each test case in the sensitivity analysis study.
As shown in Fig. 5a, while the height of the boundary surface is uniform in the cylindrical phantom, the mouse boundary is definitely nonuniform. In Fig. 5b, we see that excitation fluence peaks over the mouse target are to 20% higher than the targets in the cylinder, due to the smaller size of the mouse cross section. In Fig. 5c, we see the local geometry near the source strongly influencing the profile of the emission fluence from the fluorescent target located just below the source. The detectability from a target located below a hump is about ten times higher than a target located under a concave depression on the mouse boundary. The profile of fluence is similar for the two positions on the cylindrical surface. While an appropriately modeled cylindrical surface to represent a small animal may be useful to determine a conservative measure of specific binding of a fluorescent target, the realistic boundary model of a small animal enables the prediction of the much wider range of detectability of a target based on its anatomical location. Thus simulations on a mouse phantom can aid in the imageability of fluorescent probes to help define bounds on “design” parameters and enhance the pace of development of targeted contrast by restricting the number of actual experiments and reducing the number of animals sacrificed.
We have proposed a method for creating digital mouse phantoms for a small animal optical imaging system design from an MRI data set by performing extensive comparison of light transport with two diffusion models on a cylindrical geometry. We have identified two principle sources of variation: (a) optical properties of the biological medium—a set of noise parameters that cannot be controlled, and (b) the boundary of the phantom that may be controlled through improving the accuracy of the geometrical model. To the best of our knowledge, this is the first reported effort at standardizing imaging system performance through the assessment of sensitivity of simulated fluence to these variations in a small animal model. The method is sufficiently general and can be extended to other optical imaging applications.
However, there are several ways in which we could potentially improve the design and accuracy of the digital phantom and extend its range of performance. An important limitation of our work is that the phantom has been created from a single MRI data set through manually chosen parameters for the image processing steps. A generalized phantom should be ideally created from multiple data sets processed through semiautomatic or automatic algorithms for boundary extraction. MRI images are prone to variability in intensities arising out of field inhomogeneity, magnet strength, and varying acquisition protocols. Automatic analysis of MRI images is an inherently challenging problem. It is, therefore, reassuring that small variations in boundary do not significantly alter the computed fluence. Our digital phantom consists of only the boundary. It would be useful to extract the anatomical structures within the mouse to create a more detailed mouse phantom. In optical imaging with anatomical priors of the brain, MRI images are used as a guide to delineate white and gray matter and CSF.19 Segmenting internal organs of the mouse from MRI data sets is a problem that is significantly more complex than extracting the mouse boundary. A detailed phantom would enable a more intricate modeling of light transport through tissue layers and allow for simulating perfusion of contrast agents and pharmacokinetics through optical techniques.
The Digimouse26 mouse atlas is represented by an isotropic voxel size of for a matrix size of . However, in the application of Digimouse to simulation of bioluminescence, Dogdas have down-sampled the atlas to a matrix size of with a voxel size of . We have limited the mesh resolution of our digital phantom to primarily due to constraints imposed by the hardware. If we moved to a paradigm where we could perform adaptive meshing of the phantom, instead of driving the FE mesh with a single resolution, we could define a fine mesh close to the source and near sharp gradients in topology to mimic rapid changes in the fluence profile and meet the criteria for accurate simulation of optical transport over a wider range of tissue types with the same hardware.
Last, the phantom needs to be tested for image reconstruction in an actual imaging system for target detection through measurements of fluorescence.
The authors thank the reviewers for useful suggestions to strengthen the content of the paper, Ravi Malladi for help during manuscript revision, Srini Rajagopalan for introducing Digimouse, and Manohar Kollegal and Amey Joshi for their contribution to the FE model development.