The near-infrared window, where intrinsic optical extinction of most biological tissues is fairly low, is a favored working frame for the efficient noninvasive imaging of disease models. A high sensitivity and specificity can be obtained by applying cost-effective optical probes targeting chosen molecular events in animals.1, 2
Current planar imaging (topographic) systems allow for a sensitive high-throughput detection of fluorophore signals in living animals. They deliver a robust signal in a short acquisition time at low cost. However, they still suffer from two severe limitations: the spatial resolution and accuracy become poorer with increasing depth of the chromophore, and an absolute quantification is not at hand.3 To overcome these drawbacks, optical fluorescence tomography (OFT) systems have been proposed and implemented.4, 5 However, despite these technological developments, reports on tomographic molecular imaging in the brain of living mice remain sparse.5 Tomographic systems currently on the market require placing the animal in a matching fluid,4 which would demand the supplemental effort of ventilating the animals intended for brain studies. This justifies the striving for full noncontact technology, although the latter faces the challenge of the pigmented fur of the animal strains used to study brain diseases.
In the present paper, we propose a tomograph that is optimized to noninvasively retrieve the distribution of fluorescent markers in the brain of mice, and this special task makes it differ from the other published techniques on the technical as well as on the postprocessing level. Our system uses a fiber-based source scheme and performs noncontact detection of fluorescence photons on the dorsal side of the head. The small distances involved in the light propagation are typical for this setup and constitute a major challenge that requires an intricate acquisition scheme and an elaborate handling of the photon transport6, 7 that we solved with the Monte Carlo method.
The paper starts with the mathematical formulation of the tomographic problem (Sec. 2). In analogy to expressions developed in the diffusion regime,4, 8 a linear system of equations can be applied. Section 3 describes the setup and the specific postprocessing routines. The technique is successively tested on a computer simulation, a resin phantom, and a pigmented furry mouse.
Linear Forward Model in the Transport Regime
The setup (see Fig. 1 ) consists of light source fibers in direct contact with the scalp, which can be in the field of view of the camera. The setup thus includes source detector combinations with small separations, requiring us to treat the photon transport in the transport regime9 rather than in the diffusion regime. Starting from a hypothetical fluorochrome concentration in the tissue, we derive the governing equations of the tomographic problem.
The fluorescence radiant flux (watts) emitted from the ’th voxel at the position inside the tissue depends on the local incoming excitation flux (watts), on the local fluorochrome concentration (moles per liter), on the specific extinction coefficient (moles per liter per meter) of the fluorochrome, on the characteristic voxel size (meters) and on the quantum efficiency of the fluorochromeis dependent on the photon flux at the source position and the dampening of light caused by scattering and absorption expressed by (watts) of this fluorescence flux reaches a detector at located on the boundary of the tissue. By summing over all voxels, a nonfocal distribution can be addressed2, 3 introduce two Green’s functions , and , for the source at and the detector at . We evaluate these Green’s functions using a Monte Carlo simulation, described in detail in Sec. 3.2.1.
To adapt Eq. 3 to a specific imaging hardware, some characteristic factors such as the laser power, sensitivity of the detection system, and coupling coefficients to and from the tissue have to be taken into account. To accommodate the unknown factors, measuring the excitation flux to the detector and building the ratio have also been suggested.5, 10, 11 The factors then cancel out under the heuristic assumption that they are only slightly dependent on the wavelength. Using this ratio as input to a reconstruction process, Eq. 3 must also include a theoretical estimate of the computed tissue transmittance between the source and the detector . We can thus write. The system of uncoupled linear equations to be solved reads is called the sensitivity matrix and is an index labeling the individual source detector combinations. The calculation of is detailed in Sec. 3.2.1. Inversion of the system gives the fluorescence yield , which is proportional to the concentration of the fluorochrome in the tissue.
The tomograph sketched in Fig. 1 uses sequential light illumination by fibers located around the mouse head. Detection occurs on the top of the object via a charge-coupled device (CCD) camera. The 24 illumination fibers are inserted inside the fiber holder shown in Fig. 2a . Clamping metal sheaths act as leads and stiffeners for the flexible fibers, which are distributed in three vertical planes with adjustable interplane distance. The animal is stereotactically fastened in a custom animal slide that can be inserted sequentially into the fiber holder and in a magnetic resonance tomograph (Bruker 7-T Pharmascan; Biospin, Rheinstetten, Germany) without any motion of the animal. For optical imaging the fiber-holder and the detecting optics are placed in a dark chamber.
The illumination light is provided by laser diodes coupled to the input fibers of an optical switch (Lightech Inc., London, England). The end of the 24 output fibers ( core) are fixed to the metal clamps. The switch controls the sequence of the fiber illumination during the tomographic acquisition. We used two laser diodes at (animal phantom) and (resin phantom) to respectively excite the fluorophores Cy5.5 (Amersham Biosciences, Buckinghamshire, England) and DY751 (Mobitec GmbH, Göttingen, Germany).
The sources are disposed all around the animal head to provide an illumination that would possibly diffuse into each voxel of the head of the animal and would be minimally affected by the presence of hair. The brain of the lying animal is the targeted organ, and this fact naturally dictated the position of the detection hardware. The CCD camera (Vers’Array 512B, Roper Scientific, Trenton, New Jersey) is placed on top of the chamber, and the focusing optics are inside the chamber. A Nikkor f/1.2 manual lens focuses the upper side of the mouse head on the cooled chip of the camera, covering a roughly pixels wide region out of the pixels available. Filters are screwed to the lenses to reject the excitation light during fluorescence measurements (three interference filters, full width at half maximum ; LOT-Oriel, Darmstadt, Germany). The fluorescence light was detected at when using the DY751 fluorochrome and at when using Cy5.5 (animal phantom). To acquire the excitation images, the interference filters were replaced by neutral density filters (Melles Griot, Albuquerque, New Mexico).
Data processing prior to reconstruction
The acquired data need processing before their use in the reconstruction step. We correct for the imperfections of the system and define the detectors according to the geometry of the object of study.
Extension of the dynamic range
For images where the light source is within the field of view of the camera, a large dynamic range is required. The signal of the excitation field is orders of magnitude stronger near the fiber tip than a few millimeters away. To overcome the 16-bit limitation of the camera analog/digital converter, the excitation field images are merged from a short-illuminated image that shows an unsaturated representation of the peak region and a set of intentionally saturated images. The latter can show a blooming effect around the peak at the fiber tip but allows a good representation of the weak signal in the pixels further away. Using more than two images, the saturated parts are incrementally supplemented through the signal of images showing less blooming or less saturation after the necessary scaling to equalize the exposure time (the unsaturated image serves as the reference). This procedure, illustrated in Fig. 3 , leads to a composite image with a higher signal-to-noise ratio and a dynamic range reaching over five orders of magnitude at the price of a longer measurement time (six images were combined in the presented case).
Definition of the detectors
The light emitted from the animal head is detected in a noncontact manner using a downward-looking CCD camera that has been placed above the animal's head. Similarly to the procedure proposed by Schulz 12 we can project a regular grid on the upper surface of the object and thus define a set of regions. The pixels inside each region are binned to form an individual spatially extended detector as will be briefly described. Because the numerical simulations are a time-consuming process, we chose not to use each pixel as a detector. The number of detectors is chosen by the operator, and this defines, most of the time, a noninteger binning factor. The size of the image is expanded to the least common multiple using a cubic spline interpolant (care is taken to preserve the energy in the matrix) to allow for a classical binning procedure.
Defining detectors from the image is straightforward in the case of the trivial geometry of the resin cylinder and has been extended for the case of a pigmented mouse and is shown in Fig. 4 . The detectors (white square) can only be defined inside an operator-defined region of interest (blue polygon). Regions containing hair or being occulted by the fibers are discarded.
Especially for our measurements using the Cy5.5 dye, fluorescence images are contaminated by a fraction of the excitation field. This leakage can be corrected by subtracting from the fluorescence images a given percentage of the excitation image. This last factor is identified using a fitting method prior to the measurements.5
Calculation of sensitivity matrix
Calculation of the sensitivity matrix [matrix of Eq. 5] involves the calculation of light propagation in a geometric configuration of the medium of interest. We use a Monte Carlo transport simulation (code originally written by Barnett and further developed by Stott and Boas13 to calculate A. The simulation is carried out on a uniform three-dimensional grid raster of cubic voxels (quantitatively for the diverse experiments of this paper, see Table 1 ).
The optical properties used for the forward simulation of the resin phantom and the in vivo measurement.
|μa in 1∕mm||μs in 1∕mm||g||μa in 1∕mm||μs in 1∕mm||g|
|Resin Phantom||not applicable||0.01.||15|
The photon packets are sent from sources located on the surface of the object, and the photon flux is registered for each committed voxel of the medium using the variance reduction scheme. Scattering is handled with the Henyey-Greenstein phase function. The photon fluence in the voxels of the object for a source at , can be considered a discrete approximation of the two-point Green’s function associated with this source for each voxel, the term in Eq. 2 is its normalized form. Square detectors are placed on the upper surface of the medium. Their position and pointing direction are matched to the experiment, and their emitting angle is defined through a numerical aperture. For each source, the program delivers an approximation of the Green’s function , inside the tissue and the integrated transmittance defined in Eq. 5. Calculation of the Green’s function in Eq. 3 is performed by virtue of the reciprocity theorem.14 Thus a simulation is run in which the detectors are used as sources. The forward model involves 5 million photon packets for each of the source and detector positions.
The program allows the definition of different types of tissue. For each, we assign typical optical properties (absorption and scattering coefficient, asymmetry factor, and index of refraction). For the mouse head, we assume different optical properties for the brain and the rest of the tissue as depicted in Table 1.
Solving the inverse problem
The inversion of the ill-posed matrix can be performed by means of several algorithms10 to find the fluorochrome distribution inside the medium from measurements taken on the outside of the medium. We used the simultaneous iterative reconstruction technique (SIRT),15 which is well-documented, robust, and of straightforward in its implementation. We added a nonnegativity constraint. We used a regularization scheme to help stabilize convergence of the ill-conditioned problem.
The results of the Monte Carlo simulations easily exceeded the capacity of our computer, so we had to reduce the grid size to for all the tests we will present in this paper.
For the first test (simulated dataset, see Sec. 3.3.1), the regularization parameter values were set between 0.1 and 0.5. The reconstruction typically needed 5000 SIRT iterations (regularization factor: 0.5), which took about on a Pentium 4 personal computer, with random access memory.
Reconstructions of the experiments (of Secs. 3.3.2, 3.3.3) were performed with 50 SIRT iterations. The use of a relaxation parameter helped contain the divergence problem, but it did not show a steady behavior along the types of experiments. Ideally, the factor would be kept as close to one as possible to accelerate the convergence, but its value was bound to the amount of noise in the data to be reconstructed. A value of 0.1 for the resin phantom and of 0.001 for the mouse experiment proved to be the best trade-offs between speed and convergence.
The reconstructed distributions still proved broad (in comparison to the 1-voxel large target of the first test model, for example), so we refined them with a descent method (conjugated gradient least square method16) known to converge steadily but at a lower rate. We used 2500 iterations with a regularization factor of 0.5.
Three test models with increasing levels of difficulty were defined. First, a completely simulated dataset based on the mouse head anatomy was produced to test the setup. This corresponds to the least challenging conditions, with well-defined optical properties and a controlled noise level. The second model consisted of measurements on a resin phantom, thus, including the experimental noise and positional errors, yet still having a good estimate of the optical properties. The final evaluation was performed on augmented mouse with a small fluorescent inclusion in the brain, bringing the full complexity of future fluorescence tomography molecular imaging experiments.
Numerical model of the mouse head
A mouse head geometry was derived from magnetic resonance (MR) images (7-T Pharmascan T2 sequence, voxel size , Bruker). The brain was isolated from the rest of the tissue through a binary segmentation, and optical properties were assigned to both tissue types (see Table 1). We used 70 detectors on the top of the mouse’s head and added the 24 sources following the excitation-detection scheme presented in Fig. 1. The origin of the fluorescent signal was set by arbitrarily assigning a nonnegative fluorescence yield value to a single voxel inside the right hemisphere.
The Monte Carlo simulations ran on a medium of voxels and generated a set of simulated measurements that we corrupted with diverse amounts of noise. We modeled the fluorescence acquisition by adding the readout noise and a model standing for the photon shot noise model. Gaussian distributions having a standard deviation of one count simulate the readout noise and were uniformly added to each channel. The shot noise was simulated for each detector; each time we randomly picked a value out of a Gaussian distribution of mean being 0 and variance being the square root of the detected value. Only a percentage (which we define as the “amount of noise level”) of this picked value was added to the measurement.
The resin phantom has a cylindrical geometry ( diameter, length) with optical properties of about and . A capillary (inner diameter of BK115; Kabe Labortechnik, Nümbrecht-Elsenroth, Germany) was placed parallel to the axis at a depth of about below the surface of the cylinder and extends up to the first half of the phantom. The optical coefficients were experimentally determined at the PTB (German Institute for Standards and Technology), and the asymmetry factor of the Heynyey-Greenstein phase function is an estimation.
For the measurements presented here, the capillary was filled with a aqueous solution of the DY751 Dyomics dye (corresponding to an amount of ). The Monte Carlo calculation used a -voxel medium ( resolution), 24 sources, and an detector grid.
To simulate a fluorescent source in the frame of a brain disease molecular imaging experiment, we implanted a capsule containing of Cy5.5-N-hydroxysuccinimidyl esther (in isotonic solution) into the brain of a pigmented furry mouse (strain C57BI6). Surgery was performed according to Ref. 3, which also implies shaving the top of the skull. The MR scans of the animal (T2 scans, see Sec. 3.3.2) helped infer the position of the capsule and model the individual anatomy of the mouse for calculation of the sensitivity matrix. The optical properties are calculated according to the model of Alexandrakis 17 for and , except for the scattering coefficient of the brain, which is modeled according to the formula proposed by Bevilacqua .18 For the absorption model, the spectra of the hemoglobins were taken from Ref. 19. The model relies on an addition of elementary spectra whose weights are determined by choosing the percentage of fat, muscle, blood, bone, and skin. The composition of the mixed tissue was arbitrarily set to 80% of muscle (anisotropy factor ), 10% of bone (0.8), 5% of skin (0.8), 2.5% of fat (0.77), and 2.5% blood (0.99), leading to absorption coefficients comparable to those found in Ref. 20. The model calculation was performed for voxels ( resolution), 24 sources, and 64 detectors.
Characterization of the reconstruction quality
The quality of reconstruction was evaluated using two criteria.21 The qualitative indicators compare the position and spread of the reconstruction to the known target fluorochrome distribution. The first parameter is the object centroid error (OCE), an evaluation of the (three-dimensional) distance between the real locus of the fluorescent object and the center of gravity of voxels with a value exceeding 50% of the maximum reconstructed value. The second parameter, the spatial correlation, is an indicator of the similarity of the reconstructed to the target distribution21 and is calculated on a two-dimensional basis using(respective of ) is the distribution of the reconstructed (respective of the target) distribution over the voxels. Quantities annotated with or an upper bar are the standard deviation and the mean of the associated distributions.
First Evaluation: Numerical Simulation
Figure 5 shows the color-coded normalized fluorochrome distribution reconstructed from a simulated experiment. The fluorochrome concentration was assumed in a single voxel in the brain (white frame in the image). As expected, the reconstructed fluorochrome distribution is broader than a single voxel but the peak concentration shows a good spatial coincidence with the target voxel. Thus the OCE was measured to be lower than half a voxel .
To evaluate the spatial accuracy inside the mouse head achievable by this setup, we reiterated the preceding test by positioning the fluorochrome inclusion at 36 different voxels in a single plane. Three noise levels were evaluated as shown in Fig. 6 (0%, 0.1%, and 1% from top to bottom). The top row corresponds to a measurement without noise and thus depicts the best reachable reconstruction parameters. The OCE representing the distance in voxels between the “true” and the predicted peak location was below . The correlation coefficient, showing the sensitivity to artifacts and spreading of the reconstructed map compared with the true distribution, had its higher values in the upper part of the slice for noiseless data. Fluorescent inclusions located in the upper cortex would then be reconstructed as a volume of acceptable spatial extent. The reconstructions of inclusions in the lower part of the slice were observed to have a larger spatial extent, thus lowering the correlation coefficient. Exceptions occurred when the inclusion was close to a source. The correlation then once again showed high values due to the pointed shape of the sensitivities at these loci.
Second Evaluation: Resin Phantom
Figure 7 shows the normalized fluorochrome concentration obtained for the measurement on the resin phantom. The noise level was lower than the 0.1% case of Sec. 4.1. The body of the phantom is shown in gray and the position of the capillary containing the fluorochrome as a white trace. The reconstructed distribution is centered on the fluorescent rod, displaying its higher values at the central excitation plane. The reconstructed fluorochrome distribution is rather homogenous under the detector array (blue line), but decreases rapidly in the outer region.
Third Evaluation: Animal Phantom
Due to the much higher absorption of the animal tissue, the measurements showed a higher noise level than the preceding case. The reconstructed concentration surrounds the capsule, the OCE having a value about . The normalized fluorochrome concentration is shown in Fig. 8 . Note that the longitudinal displacement of the reconstruction is lower than in the sagittal plane, because the reconstructed concentration tends be pushed toward the boundary of the tissue.
The tomographic approach proposed here was engineered for noninvasive brain imaging in rodents. It differs from other approaches 4, 8, 12, 22, 23 by means of the simultaneous acquisition of a reflectance and transmission geometry. To enhance the spatial accuracy in the brain region, we combine the remission mode with an arrangement of sources around the head that should maximize the photon visitation inside the brain (as the detectors are located above the parietal cortex). This prerequisite demands the application of the transport equation rather than the diffusion equation and the use of arbitrary shaped geometries for the photon transport model. We will now critically discuss the advantages and limitations of this approach and propose ways to amend it.
Advantage of Tomographic Setup
The OCE of the numerical studies showed a good matching between the target and reconstructed fluorophore distribution in all three models (see Figs. 5, 7, 8). The positions of the sources and detectors were chosen for the banana-shaped sensitivities to efficiently cross the brain, albeit at the cost of a sparser presence in the lower part of the rodent head. Positioning the detectors on the apex of the head also led to the smaller spatial spread of the reconstruction on the upper part of the brain, as revealed by the higher values of the correlation coefficient in this region.
Limitations of the Proposed Geometry
The last feature also revealed a drawback. As detailed in the experimental section, the processing of the acquired images thus required an additional step for the fibers lying in the field of view of the apparatus. The ratios of Eq. 6 would not be thoroughly related to the measurements without the extension of the dynamic range of the excitation field images.
For the same reason, small distances between sources and detectors require handling the photon propagation by means of the radiative transfer theory.6, 7, 9 Using the diffusion approximation would save time for calculation of the forward model, but at the expense of spatial accuracy.
Uncertainty Sources: Optical Properties and Noise
We sent segmented MR anatomical images into the Monte Carlo algorithm, enhancing the reality level for the sensitivity matrix calculations, but we limited the number of tissues, in other words, the set of optical coefficients, due to a poor knowledge of the tissue optical properties. The values found in the literature showed a significant spread,20 reaching in some cases up to a factor 3. Despite our precaution, further studies confirmed the tendency suggested by the first test that our model was much more sensitive to noise than to errors in the absorption as in the scattering coefficient. Surprisingly low effects were observed on the localization of the target for exaggerated changes in the optical coefficients. We summarized these results in the Table 2 , based on reconstructions of a synthetic dataset similar to what was done in the first test of this paper. We observed earlier that distortions appeared in phantom reconstructions much quicker when the amount of noise was increased, pushing the reconstruction toward the border on the medium. Similarly, the numerical simulations (Sec. 4.1) and phantom measurements (Sec. 4.2) presented in this paper allowed for rather good positioned reconstructions, but a spatial offset appeared in the animal study that was inherently more prone to noise contamination. Still, for the sake of a complete tomographic scheme, one would identify the spatially resolved optical properties of the medium on the basis of preliminary measurements before trying to retrieve the fluorophore concentration.
Effect of changes in the optical properties on the reconstructed position in numerical case similar to the first test we presented.
|Case||Variation in μa||Variation in μs||OCE in mm||Spatial Correlation|
In this paper we applied our tomographic schemes in rather controlled conditions. The experimental part of the paper resulted in satisfactory reconstructions. They were nevertheless expressed in arbitrary units. We can consider these measurements as calibrations for the next experiments, as we can now relate the reconstructed value to units of concentration.
Enhancement of the technique can primarily be expected in the forward model and in the reconstruction scheme. We aim to reduce the computational effort and use more refined inversion strategies.
Extending the number of detectors means extending the computational time. As an example, the forward model of the third test was computed in about (8 virtual computer processing units). The detector distribution could possibly be optimized in the postprocessing procedure for an optimal binning scheme of the CCD images. We envision reducing computational effort and also further stabilizing the solution of the inverse problem by tailoring the detector distribution according to the information content.
As a final remark, the reconstructions were rapidly performed with the SIRT scheme, which is only a subtype of the algebraic reconstruction techniques. They exist in great variety,15 but their common semiconvergence issue led to the extensive testing of other methods 10, 11, 16, 21, 24, 25 that might perform better in our case too.
The tomographic technique presented here targets the reconstruction of fluorochrome distributions in the head of mice. The task is addressed by a special disposition of the sources and detectors in the experimental setup. The forward model of the setup is based on Monte Carlo calculations for the propagation of light in the inhomogeneous tissue and uses an anatomical input for the computation volume.
The setup was tested in three different models. The reconstructions of forward simulations showed that the system is able to meet the demands for good spatial resolution in the brain and a medium resolution in the extracerebral tissue. In a second text, using a phantom with known optical properties, the system is able to localize a fluorescent inclusion with an error of about the voxel size. In the final application, an animal measurement where the background optical properties could only be estimated, the inclusion was located with a spatial mismatch of about .
This work was supported by funding from the EFRE EU, the Bundesministerium für Bildung und Forschung, and the Hermann and Lilly Schilling Stiftung. We gratefully thank Susanne Müller for her technical support with MR imaging.