1. Introduction

The function of molecular imaging is to help study biological processes in vivo at the cellular and molecular levels [1, 2, 3, 4, 5]. It may non-invasively differentiate normal from diseased conditions. The imaging of molecular signatures, specific proteins and biological pathways allows early diagnosis and individualized therapies, leading to the so-called molecular medicine. While some classic techniques do reveal information on micro-structures of the tissues, only recently have molecular probes been developed along with associated imaging technologies that are sensitive and specific for detecting molecular targets in animals and humans. A molecular probe has a high affinity for attaching itself to a target molecule and a tagging ability with a marker molecule that can be tracked outside a living body. Among molecular imaging modalities, optical imaging, especially fluorescence and bioluminescence imaging, has attracted remarkable attention for its unique advantages, especially performance and cost-effectiveness.

Among various optical molecular imaging techniques, bioluminescence tomography (BLT) [6, 7, 8] is an emerging and promising bioluminescence imaging modality. In contrast to fluorescent imaging, there is no background auto-fluorescence with bioluminescence imaging. The introduction of BLT relative to planar bioluminescence imaging can be, in a substantial sense, compared to the development of x-ray CT based on radiography. Without BLT, bioluminescence imaging is primarily qualitative. With BLT, quantitative and localized analyses of a bioluminescent source distribution become feasible in a mouse.

Optical molecular imaging of mice is of paramount importance because of the availability of genetically homogeneous inbred strains of mice and the creation of transgenic strains carrying activated and inducible forms of oncogenes or knockouts of tumor suppressive genes. Mice are used in over 90% of mammalian research studies. Mouse models are established for over 90% of human diseases. Currently, results from leading groups all suggest that in favorable cases with strong prior knowledge BLT can produce valuable tomographic information [7, 9, 10, 8, 11, 12, 13, 14, 15, 16, 17, 18, 19]. One category of the favorable cases is that a relatively small source permissible region is known prior to the BLT reconstruction. Nevertheless, it is not always reliable or feasible to define such a permissible region effectively. In brief, it remains extremely challenging to achieve a significantly better and consistently stable BLT performance.

Zhao et al. have recently reported that bioluminescent spectra can be significantly affected by temperature [20]. As it is well known, when the luciferase is combined with the substrate luciferin, oxygen and ATP, biochemical reactions occur to emit bioluminescent photons. Luciferase enzymes from firefly (FLuc), click beetle (CBGr68, CBRed), and Renilla reniformins (hRLuc) have different emission spectra that are temperature dependent. When temperature is increased from 25°C to 39°C, the brightness increases in general. Moreover, FLuc has a 34 nm spectral red shift, as shown in Fig. 1 [20].

 

Fig. 1. Spectra of FLuc at 25°C, 29°C, 33°C and 37°C (data from [20] with permission).

Download Full Size | PDF

The Fluc spectrum may be partitioned into the following three bands: [500,590], [590,625], and [625,750] nm. The first spectral interval covers cyan, green and yellow. The second interval is essentially for orange. The last interval is for red. The rationale for this particular partition is that these intervals contain very similar amounts of bioluminescent energy at the mouse body temperature 37°C. Table 1 lists percentages of energy in each spectral band along with the normalized total energy at different temperatures. While the original data in the temperature range [25°, 39°]C are available [20, 21], the data at 41°C are from a third order polynomial extrapolation with least-square fitting.

Table 1. Percents of energy in each spectral band along with the normalized total energy at different temperatures (The data from 25°C to 39°C are from [20, 21] with permission). The data at 41°C are from our extrapolation. The total energy is normalized by the total energy at 37°C.

View Table | View all tables in this article

We have recognized the aforementioned temperature dependence of bioluminescence as a major opportunity to overcome the inherent illposedness of bioluminescence tomography (BLT). As shown in Table 1, from 25°C to 41°C the total emitted energy is increased as much as about 80%. While the energy percentages in [590,625] nm and [625,750] nm become gradually larger, that in [500,590] nm is accordingly reduced. More interestingly, when the temperature is elevated from 37°C to 41°C, the total energy in [590–750] nm will be increased by about 17%. In bioluminescence imaging, such a magnitude of a bioluminescent source change can be reliably recorded on the body surface of a mouse. For that purpose, we need to heat a bioluminescent source region in a well-controlled manner. If such a heated target is small, we will actually be able to pinpoint a source location and quantify it assuming a known physical model of the mouse. If the heated volume is not that small, we can at least delimit a much-reduced bioluminescent source permissible region, which will in turn greatly improve the BLT reconstruction quality.

In this paper, we propose to develop the first temperature-modulated bioluminescence tomography (TBT) system [22]. Our fundamental idea is to produce a small hot spot (or a differently shaped volumetric region) in a living mouse, and measure the signal of difference on the body surface of the mouse for BLT. The key components of a TBT system include a heating system such as a focused ultrasound array and a traditional BLT system. It is the temperature-modulated, spectrally and temporally resolved data collected using the TBT system that allows us to convert the ill-posed BLT problem into a better-conditioned or well-posed TBT one, and promises major gains in the BLT reconstruction performance in terms of localization, quantification and robustness. In the second section, we describe the ultrasound heating principles and a design of a focused ultrasound array. In the third section, we formulate the temperature-modulated bioluminescence imaging process and present a TBT algorithm. In the fourth section, we report our numerical simulation results. In the last section, we discuss relevant issues and conclude the paper.

2. Ultrasound heating

2.1. System design

Ultrasound has been widely used for heating in hyperthermia. In non-elastic media such as water and tissue, ultrasound energy of low intensity can be directly converted into heat by inelastic scattering. More dramatic heat production is due to the cavitation [23, 24]. Although diagnostic ultrasound normally operates at the frequency range of 1–20 MHz [25], the penetration depth of high frequency ultrasound waves is very short. Therefore, frequencies of 1–3MHz are commonly used to achieve up to 60 mm penetration depth [26] in tissue. Typically, ultrasound is generated using transducers made of piezoelectric materials. The resonance frequency of an ultrasound transducer depends on its size, shape and thickness, and can be modeled by plate and beam theory [27, 28]. Many ultrasound heating devices are constructed under the name of high intensity focused ultrasound (HIFU) [25, 29, 30, 31, 32, 33, 34, 35, 36, 37], which are made in a disk shape to focus into a small volume. The HIFU systems are mostly designed for large areas (a few centimeters in each dimension) and high temperatures (up to 80°C), see Ref. [38] for a review.

Hyperthermia has a cytotoxic effect on the cells. It causes responses in cell membranes, cytoskeleton, cellular proteins, nucleic acids and cellular immune systems [39]. After exposure to greater than 42°C for sufficiently long time, cells develop thermotolerance as an antagonist of the hyperthermic cell death. Hyperthermia also enhances the cytotoxicity of various antineoplastic agents. It induces alternations of blood flow and changes in microenvironment. An internal temperature above 50°C (122°F) will lead to immediate cell death. In this TBT project, we heat each tissue spot up to 43°C for a relatively short time period.

To improve the heating locality, a cylindrical transducer array was designed by Lu at al. [40]. This array was mounted on a 3D translation table. Numerical simulation of similar systems was reported by Ju at al. [41]. However, these applicators were intended for breast hyperthermia. Their control parameters are unsuitable for TBT of small animals. Only recently, Singh et al. designed a small animal hyperthermia ultrasound system (SAHUS) to study tumour hypoxia [42, 43, 44]. The SAHUS consisted of an acrylic applicator with up to four 5MHz ultrasound transducers, a 10-channel RF generator, a 16-channel thermocouple thermometer, and a PC. The system allows real-time temperature feedback control of power deposition. It can produce hyperthermia within a narrow temperature range (regulation error: 0.02±0.01°C) around 41.5°C. Although this system uses non-focused ultrasound transducers, it can heat small tumors of 8 mm in diameter with the aid of interstitial thermocouples. Clearly, it is desirable for TBT to produce a smaller heated volume without interstitial thermocouples.

 

Fig. 2. The cylindrical ultrasound transducer array system. The surface of each transducer is a part of a large sphere of radius R, centered at the origin of the coordinate system.

Download Full Size | PDF

As shown in Fig. 2, we propose to use a cylindrical focused ultrasound array system to demonstrate the feasibility of TBT. The ultrasound applicator is a multi-transducer array mounted on a cylindrical supporter of 10 cm diameter. The applicator may be mechanically moved with respect to the mouse so as to scan a sufficiently large volume of interest with a precision of 0.2 mm in each dimension. The spherical focusing transducers are driven by multichannel RF generators containing function generators and power amplifiers under PC control.

2.2. Pressure field expression

The ultrasound absorption power density can be estimated using an empirical formulas based on the exponential decay law of ultrasound intensity in tissue [45]. The specific absorption rate can be computed from the ultrasound pressure distribution. The ultrasound wave propagation in lossy media is usually associated with dispersion and acoustic attenuation, typically exhibiting a frequency dependency [46]. A more sophisticated treatment involves a time-domain attenuation model and multiple relaxation channels [47, 48, 49, 50, 51, 52, 53, 54, 55].

Ultrasound motion is governed by the wave equation, which gives rise to the Helmholtz equation in the frequency domain via the time harmonic approximation [56]. The full wave modeling of therapeutic ultrasound propagation in fluids was provided by Ginter et al. [57]. Alternatively, the equation for the sound pressure can be derived from the Navier-Stokes type of equations of gas dynamics. Using the conservation laws of mass and momentum, Liebler et al. [58] derived hyperbolic equations for the ultrasound pressure field. McGough et al. [59, 60] developed an efficient grid sectoring method for calculation of the near-field pressure.

In our system of N_t spherical focusing transducers targeting the center of the cylindrical ring, the pressure field at a given position x can be estimated from the Rayleigh-Sommerfeld radiation integral equation [61, 62]

where S′ denotes the total surface area of the transducers, ρ the density of the medium, c the phase velocity of the sound wave, u the velocity amplitude, λ the wavelength, k the wavenum-ber, α the attenuation coefficient, and x′ ∈ S′. The surface of each spherical focusing transducer can be partitioned into N rectangular elements. For a sufficiently fine partition, the formula by Ocheltree and Frizzell [61] can be modified to calculate the pressure field at x

where x′_n = (x′_n,y′_n,z′_n) gives the coordinates of x in the local coordinate system associated with the nth element (the origin at the center of the element, the z-axis passing through the focusing center of the transducer), L_n is the distance between x and the center of the element, h_x and h_y are the lengths of the element along the x′_n - and y′_n -axes respectively, $u_n=\sqrt{\frac{2W}{\mathit{\rho cA}}}$ is the velocity amplitude with radiation power W of a transducer, and $A=2\pi R^2\left(1-\sqrt{1-\frac{d^2}{4R^2}}\right)$ the surface area of the transducer, where d is the diameter of the transducer, and R the radius of the cylinder.

2.3. Bioheat transfer model

To reduce reflection, both the ultrasound applicator and the corresponding part of the mouse body surface are immersed in water during the heating process. Then, we can assume that the density and the speed of ultrasound in two layers are sufficiently close so that reflection and refraction at the water-tissue interface are negligible [40, 63]. Furthermore, the ultrasound power deposition can be regarded as stationary.

The ultrasound thermal transport process in phantom and in vivo is modeled using Pennes’ bioheat transfer equation [64, 40, 65]

where ρ is the tissue density, T the temperature field, T_a the normal body temperature of the mouse (also the water temperature), κ the thermal conductivity of the tissue, c_t and c_b are the specific heat of tissue and blood, respectively, ω is the blood perfusion rate, Q the absorption power density given as Q = ρq, and q is the specific absorption rate. We have [65, 41]

where α is the ultrasound absorption coefficient of tissue, c the speed of sound, and p the ultrasound pressure field. At each point, the specific absorption rate can be calculated as the sum of the contributions from all of the transducers. As heating time is relatively long and the heating power is relatively low, the temperature fluctuation can be neglected. Thus, a steady state bioheat transport equation is obtained as follows:

3. Image reconstruction

3.1. BLT formulation

Since TBT is an extension of BLT, let us first give an overview of the BLT formulation in the single and multispectral cases, using the notation in [8]. Bioluminescent photon scattering dominates over absorption in the biological tissue. Hence, this photon transport process can be modeled by the following steady-state diffusion equation [66, 8]:

where Ω denotes the region of interest, Φ(x) the photon fluence [Watts/mm²] at location x, S(x) the bioluminescent source density [Watts/mm³], μ_a (x) the absorption coefficient [mm^-1], μ_s (x) the scattering coefficient [mm^-1], g the anisotropy parameter, ∂Ω the boundary of Ω, v the unit outer normal on ∂Ω,A(x;n,n′) a function depending on the refractive indices n for Ω and n′ for the surrounding medium. Since both the ultrasound applicator and the corresponding part of the mouse body surface are immersed in water in the heating process, the surrounding medium is water with n≈n′ = 1.33. In this case A (x;n,n′) ≈ 1. In our study, we assume that the bioluminescence imaging experiment is performed in a totally dark environment. The outgoing photon density on ∂Ω is:

Then, the BLT problem is to determine a source S given measurement data Q based on Eq. (6).

Using the finite-element method to solve the BLT problem, we have [8]:

Furthermore, we can write Eq. (8) as

where $\tilde{\Phi }$ represents the measurable photon density values of the finite element nodes on the boundary ∂Ω, and Φ^* the photon density values of the internal nodes, while the source vector S is divided into S^p in a permissible region Ω_s and S ^* = 0 in the forbidden region. The permissible region is where a bioluminescent light source may exist, while the forbidden region is where no bioluminescent light source could exist. Now the equation can be reduced to [8]

where B = (M _1,1 -M _1,2 $M_{2,2}^{-1}$-$M_{1,2}^T$)^-1(F _1,1 - M _1,2-$M_{2,2}^{-1}$-F _2,1).

In a BLT experiment, the output photon flux Q(x) from a mouse is captured with a CCD camera. By Eq. (7), the photon density on the mouse body surface can be obtained from Q(x). With the measurement of $\tilde{\Phi }$ denoted as Φ^m, our problem is to recover the source S^p from Φ^m based on Eq. (10). Since the measured data are necessarily corrupted by noise, it is not optimal to invert Eq. (10) strictly. We minimize the following objective function to perform a BLT reconstruction:

where W denotes weighting matrix, η a stabilizing function, ε a regularization parameter [8]. In our study, we set η(X) = X^TX. Although the regularizaion method is widely used to overcome the ill-posedness of an inverse problem, there is not a well-accepted governing theory on the optimal selection of the stabilizing function and its regularization parameter. If the problem is more strongly regularized, the solution will be more stable but the error will be larger. In general, the regularization method must be fine-tuned in the simulation and experiments.

Now, we are ready to formulate multi-spectral BLT (MBLT) as an extension of the above single-spectral formulation. First, let us assume no knowledge on the spectral profile of an underlying bioluminescent source, such as in the case of mixed bioluminescent probes with unknown concentrations. At each wavelength of interest λ_i , we have

The corresponding objective function is

where n is the number of spectral bands, W λi a weighting matrix for spectrum λi . Then, the MBLT problem is just to find minimizers for OM(Sp ). In a more general scenario, there are several bioluminescent probes working together, these probes may have different temperature-dependent spectra. To recover the concentrations of these probes, we assume that the underlying source contain k bioluminescent probes Pj with corresponding spectral weights wPiλi , j = 1,…,k. By solving the following linear system, Xj (the power for Pj ) can be recovered.

Clearly, n ≥ k is needed for a unique solution to the linear system. Note that in the case of only one bioluminescent probe in a mouse, Eq. (12) can be enhanced as [13]

where n is the number of spectral bands, W _λi a weighting matrix for spectrum λ _i.

3.2. TBT algorithm

In the TBT process, we first record bioluminescent data on the body surface of a mouse at the body temperature, and may perform a regular BLT reconstruction. Based on all the information we collect, we can determine which regions should be heated for improvement of the BLT reconstruction. For example, we may evaluate all the elements in the permissible region via a cluster analysis, which is a well-established methodology in pattern recognition to merge finite elements of bioluminescent sources into localized groups. Then, we may use minimal spheres to cover identified clusters, and heat each of such spheres to collect bioluminescent data again. If there is a significant difference between the data before and after heating, a TBT reconstruction will be performed from the difference signal for the heated volume. If any heated volume does not produce a significant change in the bioluminescent signal, the previous BLT reconstruction in that region is considered invalid, and we need to select a different permissible region. Generally speaking, any body region R in a mouse can be heated to determine if there exists any bioluminescent source, and if so we can perform a TBT reconstruction based on the difference data measured on the mouse body surface before and after heating. Assuming that the temperature-dependent spectrum of the bioluminescent probe is known, we have the following two major TBT cases that involve (i) spectrally mixed and (ii) multi-spectral datasets, respectively.

In the case of a spectrally mixed dataset, at the reference temperature (the body temperature), we have

where $S_I^p$ and $S_O^p$ are the source vectors inside and outside the heated region, respectively. After heating, we have

where t is the normalized total energy at the elevated temperature, as shown in Table 1, w′_λi denotes the spectral weight at λ_i at the elevated temperature, and w _λi and w′_λi are available in Table 1. Now, let us compute the signal difference again,

Therefore, we have the following objective function for TBT:

where W is a weighting matrix.

In the case of a multi-spectral dataset, for each spectral band we have

After heating, we have

Then, the difference signal and the objective function become

where W _λi is a weighting matrix for spectrum λ_i .

Assuming that there are a number of temperature-dependent spectra from multiple kinds of bioluminescent probes, we generally do not know the concentrations of these probes at any particular location. In this case, we can perform a TBT reconstruction for each spectral band. After the bioluminescent power in each band is obtained, we can solve a linear equation system similar to Eq. (14) for 3D reconstruction of the concentrations of these probes.

4. Simulation results

4.1. Ultrasound heating simulation

The simulation was carried out on a workstation with a FORTRAN code. The standard second order finite difference scheme was used for discretizing the bioheat equation. A resolution of 0.1 mm/MHz was used to ensure the accuracy. In our simulation, we set the average speed of sound c = 1540 m/s, the specific heat of tissue c_t = 4186 J/kg °C and the ultrasound absorption coefficient α= 8.6 fm^1.5 Np/m [67, 68, 41, 65].

 

Fig. 3. Plots of the Q distribution on the plane z = 0 for ω= 1.5 kg m^-3 s^-1, R = 5 cm and d = 0.635 cm. (a) f = 1, N_t = 17; (b) f = 2, N_t = 29; (c) f = 3, N_t = 41.

Download Full Size | PDF

 

Fig. 4. Plots of the Q distribution on the plane z = 0 with different N_t values for f = 2 MHz, w= 1.5 kg m^-3 s^-1, R = 5 cm, and d = 1.27 cm. (a) N_t = 11, (b) N_t = 19, and (c) N_t = 29.

Download Full Size | PDF

4.1.1. Absorption power density distributions

We first examined the effect of frequency as it primarily determines the focus radius in the heating process. We selected three frequencies, f = 1,2 and 3 MHz. In all the cases, we made ω= 1.5kgm^-3s^-1,R = 5cm,d = 0.635cm. For f = 1,2, and 3 MHz, the number of transducers N_t was set to 17, 29 and 41, respectively. As shown in Fig. 3, the contour plots tightly surround the focus center, showing that the focus radius was 3, 1.5 and 1 mm for f = 1, 2, and 3 MHz, respectively. That is, with f = 2 or 3 we can obtain the desired absorption power density distribution.

Then, we fixed f = 2 MHz without loss of generality and examined the effects of other parameters. The transducer array was meant to achieve constructive-interference at the focus center and destructive-interference elsewhere. Imperfect interference led to undesired sidelobes. Since the geometry was fixed, the second most important parameter is the number of transducers N_t . Figure 4 depicts the impact of N_t on the plane z = 0. Clearly, when N_t = 11 only a small region that was very close to the focus center had a desired destructive-interference, while undesired sidelobes were scattered widely. When N_t = 19 the heating locality was much improved. When N_t = 29 an ideal absorption power density distribution was obtained. Actually, for a given frequency there is a critical number of transducers N_c to achieve a desired absorption power density distribution. As the frequency becomes higher, N_c becomes larger. For f = 1,2 and 3, N_c = 15,29 and 39, respectively. On the other hand, the maximum number of transducers is limited by Interestingly, the increment in the number of transducers did not affect the size of the central peak cross section. The diameter of the central peak cross section is about four times of the wavelength. That is, although the sound wave does not propagate in specific direction, it can be effectively focused via appropriate interference.

 

Fig. 5. Plots of the Q distribution on the plane z = 0 with different R values for f = 2 MHz, ω= 1.5 kg m^-3 s^-1, d = 0.635 cm, and N_t = 29. (a) R = 4 cm, (b) R = 5 cm, and (c) R = 6 cm.

Download Full Size | PDF

 

Fig. 6. Plots of the Q distribution on the plane z = 0 with different d values for f = 2MHz, ω= 1.5 kg m^-3 s^-1, R = 6 cm and N_t = 29. (a) d = 0.635 cm and (b) d = 1.27 cm.

Download Full Size | PDF

Figures 5 and 6 show the effects of the focus length R and the diameter d. These plots were visually similar. The only very small changes were observed in the maximum intensities with respect to different R and d values.

4.1.2. Temperature distributions

Having examined the impact of various dimensional parameters on the absorption power density distribution, we investigated how the temperature distribution was controlled by the transducer power and the blood perfusion rate. It was found that W did not change the shape of Q but affected its magnitude. Consequently, the maximum temperature T _max was a function W. Moreover, W had only a scaling effect on Q, and thus its effect on the temperature increment (T _max - T_a ) is linear (T_a = 37°C). In other words, with the other parameters fixed, the temperature increment is increased twice if the power W doubled, as shown in Fig. 7.

Finally, we investigated the impact of the blood perfusion rate ω. Figure 8 shows representative results for ω= 0.5kgm^-3s^-1 and ω= 3.5kgm^-3s^-1. It was seen that the temperature peak became sharper for a larger perfusion rate. On the other hand, ω only rescaled the steady state temperature but did not change its profile shape. In particular, the smaller ω was, the weaker its cooling effect, and the less the power W would be needed, as illustrated in Fig. 8.

In Tables 2, 3 and 4, we summarize the recommended control parameters for f = 1,2,3MHz, respectively. In these cases, T_a = 37°C, the maximum temperature was set to 43°C, and the number of transducers was minimized. Also, included in the tables are different perfusion rates, focus lengths and transducer diameters.

 

Fig. 7. Plots of the T distribution on the plane z = 0 with different W values for f = 2 MHz, ω= 1.5 kg m^-3 s^-1, R = 4 cm and d = 0.635 cm. (a) W = 0.008 w and (b) W = 0.016 w.

Download Full Size | PDF

 

Fig. 8. Plots of the T distribution on the plane z = 0 with different ω values for f = 2 MHz, R = 4 cm and d = 0.635 cm. (a) ω= 0.5 kg m^-3 s^-1, W = 0.018w and (b) ω = 3.5 kg m^-3 s^-1, W = 0.025 w.

Download Full Size | PDF

Table 2. Recommended control parameters for f = 1 MHz.

View Table | View all tables in this article

While the cylindrical array design allows an excellent absorption power density distribution on the x – y plane, we now determine how the absorption power density is distributed along the z-axis. Figure 9 shows the absorption power density at the x = 0 cross section. The plot at the y = 0 cross section is identical. Although the absorption power density distribution is symmetric on the x–y plane, there is an undesired power absorption off the focus center along the z-direction.

Table 3. Recommended control parameters for f = 2 MHz.

View Table | View all tables in this article

Table 4. Recommended control parameters for f = 3 MHz.

View Table | View all tables in this article

 

Fig. 9. Longitudinal extent of the Q distribution over the plane x = 0. (a) Surface plot and (b) a contour plot for f = 2 MHz, ω = 1.5 kg m^-3 s^-1, R = 4 cm, d = 0.635 cm and N_t = 29.

Download Full Size | PDF

In Fig. 10, temperature distributions are further analyzed along the z-axis. These quantities of interest are plotted over the y–z cross section. It is shown that the size of high temperature region (T > 40°C) is less than 1 cm. Technically, it is feasible to reduce the longitudinal extent of the heated volume using a few more focusing arrays.

 

Fig. 10. Longitudinal extent of the T distribution over the plane x = 0. (a) Surface plot and (b) a contour plot for f = 2 MHz, ω = 1.5 kg m^-3 s^-1, R = 4 cm, d = 0.635 cm, and N_t = 29.

Download Full Size | PDF

 

Fig. 11. Finite element model of the mouse chest phantom.

Download Full Size | PDF

4.2. TBT reconstruction

We tested the above-proposed TBT algorithm in the numerical simulation using a mouse chest phantom. This anatomically realistic phantom was digitally built from an MRI mouse atlas (http://www.mrpath.com/visiblemouse.html). Specifically, based on a segmented MRI image volume we obtained a finite-element model of the mouse chest as shown in Fig. 11. This phantom consisted of lungs, a heart, a liver, and muscle. The appropriate absorption μ_a and reduced scattering μ′_s coefficients were assigned to these anatomical structures as shown in Table 5, which were extracted from the literature [16, 69, 70, 71].

As mentioned before, among the four kinds of luciferase enzymes: hRLu, CBGr68, CBRed, and Fluc, only FLuc has not only an increment in the emission power but also a red-shift in the spectrum with temperature elevation. As a result, in this feasibility study we only used FLuc as the bioluminescent probe. As shown in Table 1, we have the original spectral data from 25°C to 39°C, and the extrapolated data at 41°C.

In our numerical simulation, we put respectively single and double bioluminescent sources in the mouse phantom, and performed TBT reconstructions from either spectrally mixed or multi-spectral datasets. In our tests, we partitioned the spectra over the interval [500,750] nm into the following three portions: (a) [500,590] nm, (b) [590,625] nm, and (c) [625,750] nm. We further assumed that in each portion the optical parameters would remain the same with temperature change. Since the normal temperature is 37°C, we set the elevated temperature to 41°C. From Table 1, w _a,37° = 0.2844, w _b,37° = 0.3264, w _c,37° = 0.3892, w _a,41° = 0.2340, w _b,41° = 0.3335, w _c,41° = 0.4325, and t = 1.0921. In all the cases, each source had a power of 10-picowatt (about 3 × 10⁷ photons per second) and limited within an finite-element. The camera was on for 600 seconds. The measurement data on the side surface of the phantom were generated according to the finite-element-based diffusion model Eq. (8). Because the number of emitted photons follows a Poisson distribution, we added Poisson noise to these data [72]. Then, we added 5% and 10% Gaussian noise to the Poisson data respectively to reflect the non-Poisson noise/biases in the measurement process. Figures 12 and 13 show the simulated optical data and the difference signals with one and two bioluminescent sources in the phantom, respectively.

Table 5. Optical parameters of various anatomical structures in different spectral bands.

View Table | View all tables in this article

 

Fig. 12. Optical data on the side surface of the finite-element phantom containing one bioluminescent source. (a) Original optical data, and (b) the corresponding difference data.

Download Full Size | PDF

Based on the geometrical model of the mouse, its optical parameters keyed to each anatomical region, and two surface optical measurement datasets before and after heating, the TBT algorithm can be applied to produce reliable results on the location and power of the bioluminescent source. Table 6 summarizes the representative TBT results. In each case, we reconstructed the source distribution 10 times. In all the tests, the source locations were randomly chosen while the same geometrical model and the same optical properties were used. It can be seen in Table 6 that the TBT algorithm worked very well.

 

Fig. 13. Optical data on the side surface of the finite-element phantom containing two bioluminescent sources. (a) Original optical data, (b) the corresponding difference data, (c) an oblique section showing the source locations (red), reconstructed source localtion (blue) and a heated region (green), and (d) an oblique section showing the photon density distribution. Sources 1 and 2 were inside and outside the heated volume, respectively.

Download Full Size | PDF

5. Discussions and conclusion

The simulation results on ultrasound heating have confirmed the feasibility of our focused ultrasound array design. We have analyzed the impact of a wide range of physical parameters on the absorption power density and temperature distributions, involving ultrasound frequency, transducer diameter, focal length, output power, number of transducers, spatial arrangement, blood perfusion rate, and so on. Our results have indicated that the heating locality is mainly controlled by the ultrasound frequency and spatial arrangement of the transducers. Although the use of a higher frequency leads to a smaller heated volume, it requires a larger number of transducers to suppress sidelobes. Our results have provided the guidelines for minimization of the heating locality and the number of transducers for a given frequency. Such a highly localized heating capability is instrumental for our proposed TBT.

While a cylindrical transducer array has been analyzed to produce a small hot spot in a small animal, other configurations for ultrasound heating and other kinds of heating means are possible as well and may be advantageous in some sense. For example, a spherical focusing array can be used to improve the heating locality. It is also possible to use an annular phased ultrasound array system to heat a nude mouse with coupling gel (without the use of water immersion). Such a system can achieve the whole-body access, although its locality may not be as good as the cylindrical design. On the other hand, we may design ultrasound array systems for heating a transverse slice or a linear path in a living animal. For example, a low frequency ultrasound beam can heat a cross-section of thickness four times of the ultrasound wavelength. To heat a relatively large volume, radio frequency (RF) below 20 MHz and microwave (MW) above 100 MHz are feasible options too. To achieve better resolution, high frequency MW should be used. However, the penetration length of RF and MW decreases with increment in frequency. For both RF current and MW radiation, absorbed power density decreases exponentially with the depth in tissue. Consequently, it is technically impossible to achieve both fine locality and deep penetration using either RF or MW applicators. Also, ultrasound, RF and/or MW may be combined for synergy.

Table 6. Simulated TBT results in the cases of single and double sources from spectrally mixed and multi-spectral datasets in terms of source localization and power estimation errors.

View Table | View all tables in this article

It is recognized that our heating-related computation can be improved by compensating for the anatomical heterogeneity. Indeed, sound speed varies about 3 percents inside the mouse body. Other spatially variant properties are blood profusion rate, thermal conductivity, specific heat, ultrasound absorption coefficient, reflection and scattering of ultrasound waves from acoustic impedance boundaries such as at the surfaces of bone and lungs. Since BLT is typically solved using a multi-modality approach, an independent tomographic volume of an individual mouse is assumed to be already available. Hence, a more sophisticated ultrasound focusing scheme can be developed to take into account the heterogeneous anatomy of the mouse under study, which is however beyond the scope of this paper. Nevertheless, with more accurate modeling and proper calibration, the inhomogeneous effect can be effectively incorporated, leading to better performance of the TBT system. Additionally, both ultrasound heating and ultrasound imaging may be performed using a single ultrasound device. Such an integrated system would allow more precise control of the heating process and monitoring of the temperature change, which can be improved by analysis of resultant optical measurement on the mouse body surface.

Furthermore, strictly speaking all of the physical parameters are temperature-dependent, which include mass density, sound speed, profusion rate, thermal conductivity, specific heat, and ultrasound absorption coefficient. It is possible to include thermal coupling effects in our computation. Nevertheless, these effects are very small in our case since the temperature increment is only a few degrees. Hence, we have not incorporated them in this feasibility study.

By heating a relatively small region and extracting the difference signal, the permissible region for bioluminescent source reconstruction is dramatically reduced. This helps tremendously regularize an ill-posed BLT problem into a better-conditioned or well-posed TBT framework. The original BLT algorithm highly depends on the specification of the permissible region. If such a region is incorrectly outlined or made too large, the BLT reconstruction will not be accurate and stable. Fortunately, this challenge has now been effectively addressed using the TBT approach. With ultrasound heating and difference signal detection, it becomes practical to be certain if the bioluminescent probe exists in any specified small volume. Hence, the reconstruction performance of TBT can be made significantly superior to that of traditional BLT. In all the 80 simulation runs we have performed so far, it has been found that the source localization and power estimation can be satisfactorily done with on average < 1mm and < 25% errors, respectively.

After a TBT reconstruction, the bioluminescent source distribution can be further refined based on the original bioluminescent data measured on the mouse body surface at the normal temperature (37°C). Since the TBT reconstruction may contain substantial errors due to data noise, model mismatches, system biases and so on, the computed bioluminescent data on the mouse body surface based on the TBT reconstruction may not explain the original data very well. This actually presents us a major opportunity to improve the TBT reconstruction according to the original data. For example, if we assume that the spatial support and relative variation of the bioluminescent source distribution reconstructed in the TBT are quite reliable, then we can scale the bioluminescent source distribution obtained in the TBT reconstruction so that the computed surface bioluminescent measurement fits the original data optimally. More sophisticatedly, we can statistically characterize all the sources of the imperfection in the TBT process, and perform the correction based on more quantitative constraints and more specific models.

While in this paper we have focused on 3D TBT reconstruction, it is emphasized that this idea is applicable in 2D bioluminescence imaging [3, 4, 5] and computational optical biopsy (COB) [73, 74]. The temperature-dependent spectra of the bioluminescence may also be utilized in biomedical imaging based on the bioluminescence resonance energy transfer (BRET) phenomena. We believe that our methodology and technology represent a major step forward in the field of BLT, and have a broad, important and immediate applicability in bioluminescence imaging of small animals.

In conclusion, for the first time we have proposed temperature-modulated bioluminescence tomography (TBT) and demonstrated its feasibility. Specifically, we have designed an exemplary focused ultrasound array to heat one small volume of interest (VOI) each time in a living mouse, and induce a detectable change in the optical signal on the mouse body surface. Based on this type of information, the TBT reconstruction has been simulated with excellent results. Currently, we are constructing the first TBT prototype including a focused ultrasound array. Phantom experiments and mouse studies using the TBT system will be reported in the future.