Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2898173

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Biomedical Background
- 3. Study of a Diffusion Based BLT Model
- 4. Extensions
- References

Authors

Related links

J Comput Math. Author manuscript; available in PMC 2010 July 7.

Published in final edited form as:

J Comput Math. 2008; 26(3): 324–335.

PMCID: PMC2898173

NIHMSID: NIHMS192976

Weimin Han, Department of Mathematics, University of Iowa, Iowa City, IA 52242, U.S.A. Email: ude.awoiu.htam@nahw;

See other articles in PMC that cite the published article.

Over the last couple of years molecular imaging has been rapidly developed to study physiological and pathological processes *in vivo* at the cellular and molecular levels. Among molecular imaging modalities, optical imaging stands out for its unique advantages, especially performance and cost-effectiveness. Bioluminescence tomography (BLT) is an emerging optical imaging mode with promising biomedical advantages. In this survey paper, we explain the biomedical significance of BLT, summarize theoretical results on the analysis and numerical solution of a diffusion based BLT model, and comment on a few extensions for the study of BLT.

Tomography is an important branch of imaging science and technology which targets image reconstruction from indirect measurement of an object under consideration. Among its numerous applications, tomography has been the driving force in biomedical imaging. As cornerstones of modern hospitals and clinics, x-ray computed tomography (CT), magnetic resonance imaging (MRI), nuclear and ultrasound imaging are widely applied for spatial and temporal reconstructions of anatomical and functional features, generated tremendous healthcare benefits over the past decades.

Guided by the so-called NIH Roadmap, molecular imaging has been rapidly developed to study biological processes in vivo at the cellular and molecular levels ([29], [27]). While some classic microscopic and spectroscopic techniques do reveal information on micro-structures of the tissues, only recently have molecular probes been utilized along with imaging technologies to detect and image molecular targets sensitively, specifically, and non-invasively. Among molecular imaging modalities, optical imaging is most attractive because of its unique advantages, especially performance and cost-effectiveness ([8], [30], [20]). Fluorescent and bioluminescent probes are commonly used for optical molecular imaging in preclinical studies of mice and rats as models of various human diseases, as well as to a limited extent in clinical research. In this context, fluorescence molecular tomography (FMT) ([21]) and bioluminescence tomography (BLT) ([26], [28]) are emerging as complementary optical molecular tomography modes.

Given the fast pace of the development in the BLT area and the major needs for more mathematical work, we present this survey as a reference for those mathematicians who are interested in solving cutting edge inverse problems for biomedical applications. In the following, first we explain the biomedical significance of BLT in Section 2. Then, we summarize theoretical results on the analysis and numerical solution of a diffusion-approximation based BLT model in Section 3. Finally, we discuss a few extensions of BLT in Section 4.

In the post-genomic era, great efforts are being made to associate genes to phenotypes for development of systems medicine that are predictive, preventive and personalized. An important aspect of this perspective is small animal imaging that allows *in vivo* studies at anatomical, functional, cellular and molecular levels. In molecular/cellular imaging, small animal features of interest are labeled with molecular probes ([18, 30]). 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. Optical imaging methods include florescence molecular tomography (FMT) ([21]) and bioluminescent imaging (BLI) ([22]), which are most promising because of their performance and cost-effectiveness, and already successfully used to investigate tumorigenesis, cancer metastasis, cardiac diseases, cystic fibrosis, gene therapies, drug designs and so on. Particularly, bioluminescent imaging has unique capabilities in probing molecular and cellular processes, and produces superior signal-to-noise ratios with little background autofluorescence. In the March 2005 issue of the Molecular Imaging Outlook^{1)}, Contag mentioned that BLI arose out of the frustration with sampling limitations of the standard assay techniques. Also, since the genes are duplicated with the cell division, BLI is more sensitive than other techniques such as nuclear imaging in which the radioactive signal is reduced with the cell division. Piwnica-Worms underlined in the same report that BLI could be applied to study almost all diseases in every small animal model.

Dr. Wang’s group conceptualized and developed the first bioluminescence tomography (BLT) prototype which compensates for heterogeneous scattering properties of a mouse and performs quantitative 3D reconstruction of internal sources from bioluminescent views measured on the external surface of the mouse ([26, 28, 7]). BLT has now become a rapidly developing area for optical molecular imaging. The introduction of BLT relative to planar bioluminescent imaging (BLI) can be in a substantial sense compared to the development of x-ray CT based on radiography. Without BLT, bioluminescent imaging is primarily qualitative. With BLT, quantitative and localized analyses on a bioluminescent source distribution become feasible inside a living mouse

The pre-requisites for BLT are bioluminescent probes, corresponding substrates, and subsequent signal collection. Naturally-occurring luciferases exhibit emission maxima between 480 nm and 635 nm. In principle, we may use luciferases with different spectral properties to sense various biological events. Recent results in the luciferase technology have confirmed spectrally-shifted signals from luciferases in various species and/or by mutagenesis. Among the current options, combining firefly (Photinus pyralis) luciferase (λ_{max} = 562 nm) and click beetle (Pyrophorus plagiopthalamus) (λ_{max} = 615 nm) seems attractive because they utilize the same non-toxic substrate. There are also areas for further development of bioluminescence reporters that could expand the utility of bioluminescent imaging. These include isolation of novel luciferases, mutation of known luciferases, luminescence-resonance energy transfer to red-emitting fluorescent proteins, and development of luciferase substrate analogs with different emission properties. Coincidentally, the latest development in the cooled-CCD camera technology has reached the point that allows us to detect very weak optical signals such as bioluminescent signals on the mouse body surface.

We use the symbol $\Omega \subset {\mathbb{R}}^{3}$ for the domain occupied by a biological medium under consideration. The boundary of Ω is denoted by Γ, which is assumed to be at least Lipschitz continuous. Thus, the unit outward normal vector ν exists almost everywhere (a.e.) on Γ.

Light propagation in the biological medium is described by the radiative transfer equation (RTE) ([2, 19]). Denote by ${\mathbb{S}}^{2}$ the unit sphere, and let *μ _{a}* =

$$\theta \cdot {\nabla}_{x}\varphi +{\mu}_{a}\varphi ={\mu}_{s}{\int}_{{\mathbb{S}}^{2}}k\left(\theta \cdot {\theta}^{\prime}\right)\varphi (\mathit{x},{\theta}^{\prime})d{\theta}^{\prime}+q,$$

(3.1)

where *ϕ* = *ϕ*(** x**,

$${\int}_{{\mathbb{S}}^{2}}k(\theta \cdot {\theta}^{\prime})d{\theta}^{\prime}=1.$$

In applications, Henyey–Greenstein scattering kernel function is widely used:

$${k}_{HG}\left(s\right)=\frac{1}{4\pi}\frac{1-{g}_{a}^{2}}{(1+{g}_{a}^{2}-2{g}_{{a}^{S}})3\u22152},\phantom{\rule{1em}{0ex}}-1\le s\le 1.$$

Here the parameter *g _{a}* (−1, 1) is a measure for anisotropy, with

The RTE (3.1) is to be supplemented by appropriate boundary value conditions. The forward model, namely the problem of determining the function *ϕ* from the RTE and the boundary value condition with a known light source function *q*, has been theoretically studied extensively in the literature; see, e.g., [9] for results on existence and uniqueness of solutions.

Mathematically, BLT is the source inversion problem to recover *q* from optical measurement on the domain boundary Γ, utilizing detailed knowledge on the optical properties of Ω. Note that knowledge of the individualized spatially variant optical properties is critical for BLT to work effectively.

The RTE is highly dimensional and presents a serious challenge for its accurate numerical simulations given the current level of development in computer software and hardware. However, since in the range of around 600 nm photon scattering outperforms absorption in a mouse, usually a diffusion approximation of the RTE is employed ([2, 19]). The diffusion approximation of the RTE (3.1) is the following equation:

$$-\mathrm{div}(D\nabla u)+{\mu}_{{a}^{u}}={q}_{0}\phantom{\rule{1em}{0ex}}\text{in}\phantom{\rule{thickmathspace}{0ex}}\Omega ,$$

(3.2)

where

$$u\left(x\right)=\frac{1}{4\pi}{\int}_{{\mathbb{S}}^{2}}\varphi (x,\theta )d\theta ,\phantom{\rule{1em}{0ex}}{q}_{0}\left(x\right)=\frac{1}{4\pi}{\int}_{{\mathbb{S}}^{2}}q(x,\theta )d\theta $$

are the averaged quantities for *ϕ* and *q* in all the directions. Here, $D=1\u2215\left[3({\mu}_{a}+{\mu}_{s}^{\prime})\right]$, ${\mu}_{s}^{\prime}=(1-\stackrel{\u2012}{k}){\mu}_{s}$ is the reduced scattering coefficient with

$$\stackrel{\u2012}{k}=\frac{1}{4\pi}{\int}_{{\mathbb{S}}^{2}}\theta \cdot {\theta}^{\prime}k(\theta \cdot {\theta}^{\prime})d{\theta}^{\prime},$$

which is independent of ** θ**. The equation (3.2) is to be supplemented by the boundary condition

$$u+2AD\frac{\partial u}{\partial \nu}={g}^{-}\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}\Gamma $$

(3.3)

where *g*^{−} is the incoming flux on Γ, and the differential operator */ν* denotes the outward normal derivative on Γ. The appearance of the parameter *A* in the boundary condition (3.3) is to incorporate diffuse boundary reflection arising from a refractive index mismatch between the body Ω and the surrounding medium. Discussion of the value of the parameter can be found in [5, 1]. Usually, this parameter is computed by the formula

$$A=(1+R)\u2215(1-R)$$

with a directionally varying refraction parameter

$$R=-1.4399{\eta}^{-2}+0.7099{\eta}^{-1}+0.6681+0.0636\eta $$

for some refractive index *η*. In BLT applications, the measurement is

$$g=-D\frac{\partial u}{\partial \nu}\phantom{\rule{thickmathspace}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}\Gamma \text{or part of}\phantom{\rule{thickmathspace}{0ex}}\Gamma .$$

(3.4)

The BLT problem we study is then to find a source function *q*_{0} given *g*^{−} and *g* such that (3.2), (3.3) and (3.4) are satisfied. Inverse source problems in such a pointwise formulation are the subject of numerous references. A recent reference is [10], where the objective is to identify the source function as a linear combination of monopolar and dipolar sources. We comment in passing that there is a related but different problem, the so-called diffuse optical tomography (DOT), also based on the diffusion approximation, where the aim is to find optical properties (absorption and reduced scattering coeffcients) of an object from diffuse signals generated by a controllable optical stimulation and measured on the external surface of the object. Some theoretical studies on the DOT problem are reported in [3, 2, 13, 24].

It is helpful to incorporate as much known information as possible in the problem formulation so as to reconstruct the source function more accurately. We call a subset of Ω the support of the light source if the light source function is nonzero in the subset and is zero outside the subset. In applications, usually a rough bound on the support of the light source is available. Thus, we suppose Ω_{0} Ω is a region that contains the light source support. The set Ω_{0} is known as the permissible region in the literature. It is desirable to have Ω_{0} exactly the light source support. But even if Ω_{0} is larger than the light source support, knowledge of a known Ω_{0} is still helpful in reconstruction of the light source. Accordingly, the differential equation (3.2) is written in the following more precise form:

$$-\mathrm{div}(D\nabla u)+{\mu}_{{a}^{u}}=p\chi {\Omega}_{0}\phantom{\rule{1em}{0ex}}\text{in}\phantom{\rule{thickmathspace}{0ex}}\Omega .$$

(3.5)

Here *χΩ*_{0} denotes the characteristic function of Ω_{0}, i.e., its value is 1 in Ω_{0}, and is 0 in ΩτΩ_{0}.

To avoid complicated notation, we express the BLT problem as the determination of a source function *p* in the differential equation (3.5) from two boundary conditions:

$$u+2AD\frac{\partial u}{\partial \nu}={g}_{1}\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}\Gamma ,$$

(3.6)

$$AD\frac{\partial u}{\partial \nu}={g}_{3}\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}\Gamma ,$$

(3.7)

i.e., we use the symbol *g*_{1} for *g*^{−}, *g*_{3} for −*A g*, and assume the measurement (3.4) is available on the entire boundaryΓ. We can also consider the case where the measurement is available only on a part of the boundary. Note that the influx *g*_{1} is zero in a typical BLT problem where the experiment is done in a dark environment. Combining (3.6) and (3.7) we obtain a third possible boundary condition

$$u={g}_{2}\equiv {g}_{1}-2{g}_{3}\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}\Gamma .$$

(3.8)

Only two of the three boundary conditions (3.6)–(3.8) are independent. To determine the source function *p*, we may associate one of the three boundary conditions (3.6), (3.7) or (3.8) with the differential equation (3.5) to form a boundary value problem, and choose one of the remaining boundary conditions to form the inverse problem for *p*. To be definite, in the rest of the section, we choose (3.6) as the boundary condition for the boundary value problem, and use (3.8) for the recovery of the source function *p*. In other words, we study the following problem, in pointwise form.

*Given suitably smooth functions D* > 0, *A* > 0, *μ _{a}* ≥ 0,

$$-\mathrm{div}(D\nabla u)+{\mu}_{a}u=p\chi {\Omega}_{0}\phantom{\rule{1em}{0ex}}\mathrm{in}\phantom{\rule{thickmathspace}{0ex}}\Omega ,$$

(3.9)

$$u+2AD\frac{\partial u}{\partial \nu}={g}_{1}\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}\Gamma $$

(3.10)

satisfies

$$u={g}_{2}\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}\Gamma $$

(3.11)

It is pointed out in [14] that Problem 3.1 is ill-posed: (1) in general, there are infinite many solutions; (2) when the form of the source function is specified, generally there are no solutions; and (3) the source function does not depend continuously on the data (instability). Since the BLT problem has to be solved through numerical means, lack of solution stability prevents the direct use of the pointwise formulation for practical simulations. We will study the BLT problem through minimizing the mismatch between predictions from the BVP and available measurements coupled with a regularization for stabilization.

We will use standard function spaces such as *V* = *H*^{1}(Ω), ${V}_{0}={H}_{0}^{1}\left(\Omega \right)$, *Q* = *L*^{2}(Ω_{0}), *L*^{2}(Ω), *L*^{∞} (Ω), and *L*^{∞}(Γ). For the given data, we assume *D* *L*^{∞}(Ω), *D* ≥ *D*_{0} a.e. in Ω for some constant *D*_{0} > 0; *A* *L*^{∞}(Γ), *A*_{1} ≤ *A* ≤ *A*_{2} for some constants *A*_{2} ≥ *A*_{1} > 0; and *μ _{a}*

Suppose we seek the source function *p* in a closed convex subset *Q _{ad}* of the space

$${Q}_{ad}=\{q\in Q\mid q\ge 0\mathrm{a}.\mathrm{e}.\text{in}\phantom{\rule{thickmathspace}{0ex}}{\Omega}_{0}\}.$$

We may also choose *Q _{ad}* to be the subset of non-negatively valued functions from a finite dimensional subspace of linear combinations of specified functions such as the characteristic functions of certain subsets of Ω.

For any *q* *Q*, the following weak formulation of the boundary value problem (3.9)–(3.10)

$${\int}_{\Omega}(D\nabla u\cdot \nabla v+{\mu}_{a}uv)dx+{\int}_{\Gamma}\frac{1}{2A}uvds={\int}_{{\Omega}_{0}}q\nu dx+{\int}_{\Gamma}\frac{1}{2A}{g}_{1}vds\phantom{\rule{1em}{0ex}}\forall v\in V.$$

(3.12)

has a unique solution *u* = *u*(*q*) *V* by an application of the Lax-Milgram Lemma ([4, 12]). Following the idea of Tikhonov regularization (e.g. [25, 11]), we let

$${J}_{\epsilon}\left(q\right)=\frac{1}{2}{\mid \mid u\left(q\right)-{g}_{2}\mid \mid}_{{L}^{2}\left(\Gamma \right)}^{2}+\frac{\epsilon}{2}{\mid \mid q\mid \mid}_{Q}^{2},\phantom{\rule{1em}{0ex}}\epsilon \ge 0$$

and introduce the following problem which is similar to the one studied in [14].

*Find* p_{ε} *Q _{ad}*

We have the following results concerning Problem 3.2.

- For any
*ε*> 0, Problem 3.2 has a unique solution*p*_{ε}*Q*. Moreover, the solution_{ad}*p*_{ε}*Q*is characterized by a variational inequality_{ad}When$${(u\left({p}_{\epsilon}\right)-{g}_{2},u\left(q\right)-u\left({p}_{\epsilon}\right))}_{{L}^{2}\left(\Gamma \right)}+\epsilon ({p}_{\epsilon},q-{p}_{\epsilon})Q\ge 0\phantom{\rule{1em}{0ex}}{\forall}_{q}\in {Q}_{ad}.$$*Q*_{ad}*Q*is a subspace, the variational inequality is reduced to a variational equation$${(u\left({p}_{\epsilon}\right)-{g}_{2},u\left(q\right)-u\left(0\right))}_{{L}^{2}\left(\Gamma \right)}+\epsilon ({p}_{\epsilon},q)Q=0\phantom{\rule{1em}{0ex}}\forall q\in {Q}_{ad}.$$ - The solution
*p*of Problem 3.2 depends continuously on all the data._{ε} - Assume the solution set
*S*_{0}for Problem 3.2 with*ε*= 0 is nonempty (this is valid if e.g.,*Q*is bounded). Then it is closed and convex. Moreover,_{ad}where$${p}_{\epsilon}\to {p}_{0}\phantom{\rule{thickmathspace}{0ex}}\text{in}\phantom{\rule{thickmathspace}{0ex}}Q,\text{as}\phantom{\rule{thickmathspace}{0ex}}\epsilon \to 0,$$*p*_{0}*S*_{0}is the unique element with minimal*Q*-norm among the solutions of Problem 3.2 for*ε*= 0:$$\mid \mid {p}_{0}\mid \mid Q=\underset{q\in {S}_{0}}{\mathrm{inf}}\mid \mid q\mid \mid Q.$$ - If
*S*_{0}= {*p*}, then we have the convergence$${p}_{\epsilon}\to p\phantom{\rule{thickmathspace}{0ex}}\text{in}\phantom{\rule{thickmathspace}{0ex}}Q,\phantom{\rule{thickmathspace}{0ex}}\text{as}\phantom{\rule{thickmathspace}{0ex}}\epsilon \to 0.$$

For a numerical approximation of Problem 3.2, we use the finite element method to solve the boundary value problem (3.12). Let ${\left\{{\mathcal{T}}_{h}\right\}}_{h}$ (*h*: meshsize) be a regular family of finite element partitions of $\stackrel{\u2012}{\Omega}$ such that each element at the boundary Γ has at most one non-straight face (or at most one curved side when we consider a two-dimensional analogue of the BLT problem). For each triangulation ${\mathcal{T}}_{h}$, let *V ^{h} H*

$${\int}_{\Omega}(D{\nabla}_{u}^{h}\cdot {\nabla}_{v}^{h}+{\mu}_{a}{u}^{h}{v}^{h})dx+{\int}_{\Gamma}\frac{1}{2A}{u}^{h}{v}^{h}ds={\int}_{{\Omega}_{0}}q{v}^{h}dx+{\int}_{\Gamma}\frac{1}{2A}{g}_{1}{v}^{h}ds\phantom{\rule{1em}{0ex}}\forall {v}^{h}\in {V}^{h}.$$

(3.13)

Corresponding to the functional *J _{ε}* (·), let

$${J}_{\epsilon}^{h}\left(q\right)=\frac{1}{2}{\parallel {u}^{h}\left(q\right)-{g}_{2}\parallel}_{{L}^{2}\left(\Gamma \right)}^{2}+\frac{\epsilon}{2}{\parallel q\parallel}_{Q}^{2},\phantom{\rule{1em}{0ex}}\epsilon \ge 0.$$

The admissible source function set *Q _{ad}* may or may not need to be discretized. In general, let ${\stackrel{~}{Q}}_{ad}\subset {Q}_{ad}$ be non-empty, closed and convex. Later in the section, we will consider two possible choices of ${\stackrel{~}{Q}}_{ad}$. We then introduce the following discretization of Problem 3.2.

*Find* ${p}_{\epsilon}^{h}\in {\stackrel{~}{Q}}_{ad}$ *such that* ${J}_{\epsilon}^{h}\left({p}_{\epsilon}^{h}\right)=\mathrm{inf}\left\{{J}_{\epsilon}^{h}\left(q\right):q\in {\stackrel{~}{Q}}_{ad}\right\}$.

Problem 3.3 has properties similar to those listed above for Problem 3.2.

For error estimation, we assume additionally that Γ *C*^{1,1}, *D* *C*^{0,1}(Ω‾), *A*^{−1} *H*^{1/2}(Γ), and g_{1} *H*^{1/2}(Γ). We say the admissible set *Q _{ad}* and the boundary data

We distinguish two cases regarding the choice of the set ${\stackrel{~}{Q}}_{ad}$. First, with the choice ${\stackrel{~}{Q}}_{ad}={Q}_{ad}$, it can be shown that for some constant *c* > 0 independent of *ε* and *h*,

$${\parallel u\left({p}_{\epsilon}\right)-{u}^{h}\left({p}_{\epsilon}^{h}\right)\parallel}_{{L}^{2}\left(\Gamma \right)}+{\epsilon}^{1\u22152}{\parallel {p}_{\epsilon}-{p}_{\epsilon}^{h}\parallel}_{Q}\le c{h}^{3\u22154}{\parallel u\left({p}_{\epsilon}\right)-{g}_{2}\parallel}_{{L}^{2}\left(\Gamma \right)}^{1\u22152}{\parallel {p}_{\epsilon}-{p}_{\epsilon}^{h}\parallel}_{Q}^{1\u22152}+c{h}^{3\u22152}[{\parallel {p}_{\epsilon}\parallel}_{Q}+{\parallel {g}_{1}\parallel}_{{H}^{1\u22152}\left(\Gamma \right)}].$$

Consequently, if *Q _{ad}* is a bounded set in

$${\parallel u\left({p}_{\epsilon}\right)-{u}^{h}\left({p}_{\epsilon}^{h}\right)\parallel}_{{L}^{2}\left(\Gamma \right)}+{\epsilon}^{1\u22152}{\parallel {p}_{\epsilon}-{p}_{\epsilon}^{h}\parallel}_{Q}\le c{h}^{3\u22154}$$

And if *Q _{ad}*,

$${\parallel u\left({p}_{\epsilon}\right)-{u}^{h}\left({p}_{\epsilon}^{h}\right)\parallel}_{{L}^{2}\left(\Gamma \right)}+{\epsilon}^{1\u22152}{\parallel {p}_{\epsilon}-{p}_{\epsilon}^{h}\parallel}_{Q}\le c{h}^{3\u22152}$$

Next, consider the case where ${\stackrel{~}{Q}}_{ad}$ is constructed with a discretization of the set *Q _{ad}*. In addition to the regular family of finite element partitions ${\left\{{\mathcal{T}}_{h}\right\}}_{h}$ of Ω‾, let ${\left\{{\mathcal{T}}_{0,H}\right\}}_{H}$ be a regular family of finite element partitions of Ω‾

$${\parallel u\left({p}_{\epsilon}\right)-{u}^{h}\left({p}_{\epsilon}^{h,H}\right)\parallel}_{{L}^{2}}\left(\Gamma \right)+{\epsilon}^{1\u22152}{\parallel {p}_{\epsilon}-{p}_{\epsilon}^{h,H}\parallel}_{Q}\le c{\parallel u\left({p}_{\epsilon}\right)-{g}_{2}\parallel}_{{L}^{2}\left(\Gamma \right)}^{1\u22152}[{H}^{1\u22152}{E}^{H}\left({{p}_{\epsilon}}^{1\u22152}\right)+{h}^{3\u22154}{\parallel {p}_{\epsilon}-{p}_{\epsilon}^{h,H}\parallel}_{Q}^{1\u22152}]+cH{E}^{H}\left({p}_{\epsilon}\right)+c{h}^{3\u22152}[{\parallel {p}_{\epsilon}\parallel}_{Q}+{\parallel {g}_{1}\parallel}_{{H}^{1\u22152}\left(\Gamma \right)}]$$

Consequently, if *Q _{ad}* is bounded in

$${\parallel u\left({p}_{\epsilon}\right)-{u}^{h}\left({p}_{\epsilon}^{h,H}\right)\parallel}_{{L}^{2}\left(\Gamma \right)}+{\epsilon}^{1\u22152}{\parallel {p}_{\epsilon}-{p}_{\epsilon}^{h,H}\parallel}_{Q}\le c[{H}^{1\u22152}{E}^{H}{\left({p}_{\epsilon}\right)}^{1\u22152}+{h}^{3\u22154}].$$

And when *Q _{ad}*,

$${\parallel u\left({p}_{\epsilon}\right)-{u}^{h}\left({p}_{\epsilon}^{h,H}\right)\parallel}_{{L}^{2}\left(\Gamma \right)}+{\epsilon}^{1\u22152}{\parallel {p}_{\epsilon}-{p}_{\epsilon}^{h,H}\parallel}_{Q}\le c[{h}^{3\u22152}+{H}^{1\u22152}{\epsilon}^{1\u22154}{E}^{H}{\left({p}_{\epsilon}\right)}^{1\u22152}+H{E}^{H}\left({p}_{\epsilon}\right)].$$

Numerical examples can be found in [14] showing the performance of the proposed numerical methods.

In this section, we point out a few extensions of the BLT model studied in Section 3.

First, recall that the goal of BLT is to produce a quantitative reconstruction of a bioluminescent source distribution within a living mouse from bioluminescent signals measured on the body surface of the mouse. While in most BLT studies so far the optical parameters of the key anatomical regions are assumed known from the literature or diffuse optical tomography (DOT), these parameters cannot be very accurate in general. In [16], we propose and study a new BLT approach that optimizes optical parameters when an underlying bioluminescent source distribution is reconstructed to match the measured data. We prove the solution existence and the convergence of numerical methods. Also, we present numerical results to illustrate the utility of our approach and evaluate its performance.

Second, a two regularization parameter framework for the BLT problem is introduced and analyzed in [6]. Similar to the discussion in [6], for any *q* *Q*, we denote by *u*_{1} = *u*_{1}(*q*) *V* the solution of the problem (3.12), and denote by *u*_{2} = *u*_{2}(*q*) *g*_{2} + *V*_{0} the solution of the problem

$${\Lambda}_{i}=[{\lambda}_{i-1},{\lambda}_{i}),1\le i\le {i}_{0}-1,\phantom{\rule{1em}{0ex}}{\Lambda}_{i0}=[{\lambda}_{i0-1},{\lambda}_{i0}].$$

This is a weak formulation of the boundary value problem defined by (3.9) and (3.11). For fixed constants *r*_{1}, *r*_{2} ≥ 0 with *r*_{1} + *r*_{2} = 1, we define the functional

$$-\mathrm{div}\left({D}_{i}{\nabla}_{{u}_{ij}}\right)+{\mu}_{a,i}{u}_{ij}={p}_{ij}\chi {\Omega}_{j}\phantom{\rule{1em}{0ex}}\text{in}\phantom{\rule{thickmathspace}{0ex}}\Omega .$$

(4.1)

and introduce the following problem.

*Find p*_{ε,r1,r2} *Q*_{ad }*such that* *J*_{ε,r1,r2} (*p*_{ε,r1,r2}) = inf{*J*_{ε,r1,r2} (q) : q Q_{ad}}.

All the theoretical results presented in the previous section can be extended to the analysis of Problem 4.1 and its numerical approximations. Note that when *r*_{2} = 0, Problem 4.1 reduces to Problem 3.2. Numerical results reported in [6] suggest that it is beneficial to choose the two regularization parameters *r*_{1}, *r*_{2} and the finite element mesh-size *h* such that *r*_{2} = *O*(*r*_{1}*h*).

Third, let us discuss at some length a general mathematical theory for the study of multi-spectral BLT. With simultaneous use of multiple optical reporters it becomes feasible to capture and decompose composite molecular and cellular signatures under *in vivo* conditions. That is, multispectral data can be measured in spectral bands on the body surface of a mouse, and the distributions of multiple biomarkers can be reconstructed in an integrated fashion using a sophisticated algorithm. In [17], a comprehensive mathematical framework for multispectral BLT is introduced and analyzed for the most general situation of using multiple bioluminescent reporters whose spectral characteristics may be a ected by their *in vivo* environment. In multi-spectral BLT, the spectrum is divided into certain numbers of bands, say *i*_{0} bands Λ_{1}, … Λ_{i0}, with

$${u}_{ij}+2A{D}_{i}\frac{\partial {u}_{ij}}{\partial \nu}=0\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}\Gamma .$$

(4.2)

Here, λ_{0} < λ_{1} < … < λ_{i0} is a partition of the spectrum range. Let there be *j*_{0} biomarkers with bioluminescent source distributions *p _{j}χΩ_{i}*, 1 ≤

$${\stackrel{~}{f}}_{i}=-{D}_{i}\frac{\partial}{\partial \nu}\sum _{j=1}^{j0}{u}_{ij}\left({q}_{ij}\right)=\frac{1}{2A}\sum _{j=1}^{j0}{u}_{ij}\left({q}_{ij}\right)\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}{\Gamma}_{i},\phantom{\rule{1em}{0ex}}1\le i\le {i}_{0}.$$

(4.3)

Here, ${D}_{i}\left(x\right)=1\u2215\left[3({\mu}_{a,i}\left(x\right)+{\mu}_{s,i}^{\prime}\left(x\right))\right]$, *μ′ _{a,i}*(

$${\mathfrak{q}}_{\ast j}={({q}_{1j},\cdots ,{q}_{i0j})}^{T},\phantom{\rule{1em}{0ex}}{\mathfrak{q}}_{i\ast}=({q}_{i1},\cdots ,{q}_{ij0})$$

With the emission filters of bandpasses Λ_{i}, the measured quantities are the outgoing flux densities ([23]):

$$\begin{array}{c}S\left({\mathfrak{q}}_{\ast j}\right)=\sum _{i}{q}_{ij},\phantom{\rule{thickmathspace}{0ex}}{\mathcal{l}}_{i}\left({\mathfrak{q}}_{\ast j}\right)={q}_{ij}-{w}_{ij}S\left({\mathfrak{q}}_{\ast j}\right),\phantom{\rule{1em}{0ex}}\mathcal{l}\left(\mathfrak{q}\right)=\left({\mathcal{l}}_{i}\left({\mathfrak{q}}_{\ast j}\right)\right),\hfill \\ {U}_{i}\left({\mathfrak{q}}_{i\ast}\right)=\sum _{j}{u}_{ij}\left({q}_{ij}\right),\phantom{\rule{1em}{0ex}}U\left(\mathfrak{q}\right)=\left({U}_{i}\left({\mathfrak{q}}_{i\ast}\right)\right).\hfill \end{array}$$

We assume that Γ_{i} is a non-trivial part of the boundary, i.e., meas (Γ_{i}) > 0. Thus, we allow the situation where the measurement of the outgoing flux densities is available only on parts of the boundary Γ.

Let us introduce some notations to simplify the exposition. The range of the index *i* is {1, … , *i*_{0}}, and that of *j* is {1, … *j*_{0}}; in particular, ∑_{i} stands for ${\Sigma}_{i=1}^{{i}_{0}}$ and ∑_{j} stands for ${\Sigma}_{j=1}^{{j}_{0}}$. Matrix $\left({\mathbb{R}}^{{i}_{0}\times {j}_{0}}\right)$ valued variables, as well as their row or column vectors, will be indicated by Euler Fraktur alphabets, e.g., $\mathfrak{p}=\left({p}_{ij}\right)$,$\mathfrak{q}=\left({q}_{ij}\right)$, $\mathfrak{u}=\left({u}_{ij}\right)$, and

$${\stackrel{~}{f}}_{i}=-{D}_{i}\frac{\partial {U}_{i}\left({\mathfrak{q}}_{i\ast}\right)}{\partial \nu}=\frac{1}{2A}{U}_{i}\left({\mathfrak{q}}_{i\ast}\right)\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}{\Gamma}_{i}.$$

Vector valued variables are indicated by boldface math fonts. We denote

$$\mathfrak{Q}=\{\mathfrak{q}=\left({q}_{ij}\right):{q}_{ij}\in {Q}_{j}\}$$

Then the boundary measurement equation (4.3) can be written as

$${(\mathfrak{p},\mathfrak{q})}_{\mathfrak{Q}}=\sum _{i,j}{w}_{ij}({p}_{ij},{q}_{ij}){Q}_{j},\phantom{\rule{1em}{0ex}}{\parallel \mathfrak{q}\parallel}_{\mathfrak{Q}}={(\mathfrak{q},\mathfrak{q})}_{\mathfrak{Q}}^{1\u22152}$$

For a vector valued variable with a subscript, we use “_{,j}” to indicate its *j*th component, e.g., *p _{ε}* = (

Let *Q*_{j} = *L*^{2}(Ω_{j}), *G*_{i} = *L*^{2}(Γ_{i}). Denote by *Q*_{ad,j} the admissible set for *p _{ij}*. We assume

$${\mathfrak{Q}}_{ad}=\{\mathfrak{q}\in \mathfrak{Q}:{q}_{ij}\in {Q}_{ad,j}\}.$$

with the inner product and norm:

$${(\mathfrak{p},\mathfrak{q})}_{\mathfrak{Q}}=\sum _{i,j}{w}_{ij}({p}_{ij},{q}_{ij}){Q}_{j},\phantom{\rule{1em}{0ex}}{\parallel \mathfrak{q}\parallel}_{\mathfrak{Q}}={(\mathfrak{q},\mathfrak{q})}_{\mathfrak{Q}}^{1\u22152}$$

for some positive weighting constants *w*_{ij}. We seek the unknown source field $\mathfrak{p}=\left({p}_{ij}\right)$ of the multispectral BLT problem in

$${\mathfrak{Q}}_{ad}=\{\mathfrak{q}\in \mathfrak{Q}:{q}_{ij}\in {Q}_{ad,j}\}.$$

With possibly different positive weighting constants ${\mathcal{w}}_{\mathcal{l},ij}$, we let

$${(\mathcal{l}\left(\mathfrak{p}\right),\mathcal{l}\left(\mathfrak{q}\right))}_{\mathfrak{Q}l}=\sum _{i,j}{\omega}_{l,ij}{\left({\mathcal{l}}_{i}({\mathfrak{p}}_{\ast j},{\mathcal{l}}_{i}){\left(\mathfrak{q}\right)}_{\ast j}\right)}_{{Q}_{j}},\phantom{\rule{1em}{0ex}}{\mid \mathcal{l}\left(\mathfrak{q}\right)\mid}_{\mathfrak{Q}l}={(\mathcal{l}\left(\mathfrak{q}\right),\mathcal{l}\left(\mathfrak{q}\right))}_{\mathfrak{Q}l}^{1\u22152}.$$

We also need the space *G* = *G*_{1} × *G*_{2} × … × *G*_{i0}, endowed with the inner product and norm

$$(f,g)G=\sum _{i}{w}_{i}({f}_{i},{g}_{i}){G}_{i},\phantom{\rule{1em}{0ex}}\parallel g\parallel G={(g,g)}_{G}^{1\u22152}$$

with positive constants *w _{i}*.

We assume $\Omega \subset {\mathbb{R}}^{d}(d\le 3)$ is a non-empty, open, bounded set with a Lipschitz boundary Γ, *A*(** x**) [

For any *q* *Q _{j}*, define

$$\sum _{\Omega}[{D}_{i}\nabla {u}_{ij}\left(q\right).\nabla \nu +{\mu}_{a,i}{u}_{ij}\left(q\right)\nu ]dx+{\int}_{\Gamma}\frac{1}{2A}{u}_{ij}\left(q\right)\nu ds={\int}_{{\Omega}_{j}}q\nu dx\phantom{\rule{1em}{0ex}}\forall \nu \in V.$$

(4.4)

Write $f=2A{\stackrel{~}{f}}_{i}$, ** f** = (

$${J}_{\epsilon M}\left(\mathfrak{q}\right)=\frac{1}{2}[{\parallel U\left(\mathfrak{q}\right)-f\parallel}_{G}^{2}+\epsilon {\parallel \mathfrak{q}\parallel}_{\mathfrak{Q}}^{2}+M{\mid \mathcal{l}\left(\mathfrak{q}\right)\mid}_{\mathfrak{Ql}}^{2}].$$

We then introduce the following problem.

*Find* ${\mathfrak{p}}_{\epsilon M}\in {\mathfrak{Q}}_{ad}$ *such that* ${J}_{\epsilon M}\left({\mathfrak{p}}_{\epsilon M}\right)=\mathrm{inf}\{{J}_{\epsilon M}\left(\mathfrak{q}\right):\mathfrak{q}\in {\mathfrak{Q}}_{ad}\}$.

We have the following results for the problem.

- Problem 4.2 with
*ε*> 0 has a unique solution ${\mathfrak{p}}_{\epsilon M}\in {\mathfrak{Q}}_{ad}$, and the solution ${\mathfrak{p}}_{\epsilon M}\in {\mathfrak{Q}}_{ad}$ is characterized by a variational inequalityWhen$${(U\left({\mathfrak{p}}_{\epsilon M}\right)-f,U(\mathfrak{q}-{\mathfrak{p}}_{\epsilon M}))}_{G}+\epsilon {({\mathfrak{p}}_{\epsilon M},\mathfrak{q}-{\mathfrak{p}}_{\epsilon M})}_{\mathfrak{Q}}+M{(\mathcal{l}\left({\mathfrak{p}}_{\epsilon M}\right),\mathcal{l}(\mathfrak{q}-{\mathfrak{p}}_{\epsilon M}))}_{\mathfrak{Ql}}\ge 0\phantom{\rule{1em}{0ex}}\forall \mathfrak{q}\in {\mathfrak{Q}}_{ad}.$$*Q*_{ad,j}*Q*are subspaces, the inequality is reduced to a variational equation_{j}$${(U\left({\mathfrak{p}}_{\epsilon M}\right)-f,U\left(\mathfrak{q}\right))}_{G}+\epsilon {({\mathfrak{p}}_{\epsilon M},\mathfrak{q})}_{\mathfrak{Q}}+M{(\mathcal{l}\left({\mathfrak{p}}_{\epsilon M}\right),\mathcal{l}\left(\mathfrak{q}\right))}_{\mathfrak{Ql}}=0\phantom{\rule{1em}{0ex}}\mathfrak{q}\in {\mathfrak{Q}}_{ad}.$$ - The solution ${\mathfrak{p}}_{\epsilon M}$ of Problem 3.2 depends continuously on the data.
- Suppose 0
*Q*. Then as_{ad,j}*M*→ ∞, ${\mathfrak{p}}_{\epsilon M}\to {\mathfrak{p}}_{\epsilon}=\left({\omega}_{ij}{p}_{\epsilon ,j}\right)$ in $\mathfrak{Q}$, where*p*_{ε}= (*p*) with ${p}_{\epsilon ,j}=S\left({\mathfrak{p}}_{\epsilon ,\ast j}\right)$, and_{ε,j}*p*_{ε}*Q*is the unique solution of the problem_{ad}where$${J}_{\epsilon}\left({p}_{\epsilon}\right)=\mathrm{inf}\{{J}_{\epsilon}\left(q\right):q\in {Q}_{ad}\},\phantom{\rule{1em}{0ex}}{J}_{\epsilon}\left(q\right)=\frac{1}{2}[{\parallel \mathit{W}\left(q\right)-f\parallel}_{G}^{2}+\epsilon {\parallel q\parallel}_{Q}^{2}],$$=*Q**Q*_{1}× … ×*Q*_{j0}and*Q*_{ad}=*Q*_{1,ad}× … ×*Q*_{j0,ad}. - Assume ${\mathfrak{S}}_{0M}$, the solution set of Problem 4.2 with
*ε*= 0, is nonempty. Then ${\mathfrak{S}}_{0M}$ is closed and convex. Moreover,where ${\mathfrak{p}}_{0M}\in {\mathfrak{S}}_{0M}$ satisfies$${\mathfrak{p}}_{\epsilon M}\to {\mathfrak{p}}_{0M}\phantom{\rule{thickmathspace}{0ex}}\text{in}\phantom{\rule{thickmathspace}{0ex}}\mathfrak{Q},\phantom{\rule{thickmathspace}{0ex}}\text{as}\phantom{\rule{thickmathspace}{0ex}}\epsilon \to 0,$$In particular, if the solution set ${\mathfrak{S}}_{0M}=\left\{{\mathfrak{p}}_{M}\right\}$ is a singleton. Then$${\parallel {\mathfrak{p}}_{0M}\parallel}_{\mathfrak{Q}}=\mathrm{inf}\{{\parallel \mathfrak{q}\parallel}_{\mathfrak{Q}}:\mathfrak{q}\in {\mathfrak{S}}_{0M}\}.$$$${\mathfrak{p}}_{\epsilon M}\to {\mathfrak{p}}_{m}\phantom{\rule{thickmathspace}{0ex}}\text{in}\phantom{\rule{thickmathspace}{0ex}}\mathfrak{Q},\phantom{\rule{thickmathspace}{0ex}}\text{as}\phantom{\rule{thickmathspace}{0ex}}\epsilon \to 0.$$

One particular multispectral BLT problem is discussed in [15], where some numerical results are reported.

Finally, we remark that the RTE based BLT problem is being under investigation.

^{1)}Supported by NIH grant EB001685 and Mathematical and Physical Sciences Funding Program fund of the University of Iowa

^{1)}http://www.diagnosticimaging.com/molecularimagingoutlook/2005mar/02.jhtml

Dedicated to Professor CUI Junzhi on the occasion of his 70th birthday

*Mathematics subject classification*: 65N21, 92C55.

Weimin Han, Department of Mathematics, University of Iowa, Iowa City, IA 52242, U.S.A. Email: ude.awoiu.htam@nahw.

Ge Wang, Division of Biomedical Imaging, Virginia Tech–Wake Forest University School of Biomedical Engineering and Sciences, Virginia Polytechnic Institute and State University, Blacksburg, Virginia 24061, U.S.A. Email: gro.eeei@gnaw-eg..

[1] Aronson R. In: Chance B, Alfano RR, editors. Extrapolation distance for diffusion of light; Photon Migration and Imaging in Random Media and Tissues (Proc. SPIE 1888); 1993.pp. 297–305.

[2] Arridge SR. Optical tomography in medical imaging. Inverse Problems. 1999;15:R41–R93.

[3] Arridge SR, Lionheart WRB. Nonuniqueness in diffusion-based optical tomography. Optics Letters. 1998;23:882–884. [PubMed]

[4] Atkinson K, Han W. Theoretical Numerical Analysis: A Functional Analysis Framework. second edition. Vol. 39. Springer-Verlag; New York: 2005. Texts in Applied Mathematics.

[5] Case MC, Zweifel PF. Linear Transport Theory. Addison-Wesley; NY: 1967.

[6] Cheng X-L, Gong R-F, Han W. A generalized mathematical framework for bioluminescence tomography. Computer Methods in Applied Mechanics and Engineering. 2008;197:524–535.

[7] Cong W-X, Wang G, Kumar D, Liu Y, Jiang M, Wang L-H, Hoffman EA, McLennan G, McCray PB, Zabner J, Cong A. A practical reconstruction method for bioluminescence tomography. Optics Express. 2005;13:6756–6771. [PubMed]

[8] Contag CH, Ross BD. It’s not just about anatomy: *In vivo* bioluminescence imaging as an eyepiece into biology. Journal of Magnetic Resonance Imaging. 2002;16:378–387. [PubMed]

[9] Dautray R, Lions J-L. Mathematical Analysis and Numerical Methods for Science and Technology. Volume 6, Evolution Problems II. Springer; Berlin: 2000.

[10] El Badia A. Inverse source problem in an anisotropic medium by boundary measurements. Inverse Problems. 2005;21:1487–1506.

[11] Engl HW, Hanke M, Neubauer A. Regularization of Inverse Problems. Kluwer Academic Publishers; Dordrecht: 1996.

[12] Evans LC. Partial Differential Equations. American Mathematical Society; 1998.

[13] Gibson AP, Hebden JC, Arridge SR. Recent advances in diffuse optical imaging. Phys. Med. Biol. 2005;50:R1–R43. [PubMed]

[14] Han W, Cong W-X, Wang G. Mathematical theory and numerical analysis of bioluminescence tomography. Inverse Problems. 2006;22:1659–1675.

[15] Han W, Cong W-X, Wang G. Mathematical study and numerical simulation of multispectral bioluminescence tomography. International Journal of Biomedical Imaging. 2006;2006 Article ID 54390, 10 pages, doi:10.1155/IJBI/2006/54390. [PMC free article] [PubMed]

[16] Han W, Kazmi K, Cong W-X, Wang G. Bioluminescence tomography with optimized optical parameters. Inverse Problems. 2007;23:1215–1228.

[17] Han W, Wang G. Theoretical and numerical analysis on multispectral bioluminescence tomography. IMA Journal of Applied Mathematics. 2007;72:67–85.

[18] Herschman HR. Molecular imaging: Looking at problems, seeing solutions. Science. 2003;302:605–608. [PubMed]

[19] Natterer F, WÜbbeling F. Mathematical Methods in Image Reconstruction. SIAM; Philadelphia: 2001.

[20] Ntziachristos V, Ripoll J, Wang L-V, Weissleder R. Looking and listening to light: the evolution of whole-boday photonic imaging. Nature Biotechnology. 2005;23:313–320. [PubMed]

[21] Ntziachristos V, Tung CH, et al. Fluorescence molecular tomography resolves protease activity in vivo. Nature Medicine. 2002;8:757–761. [PubMed]

[22] Rice BW, Cable MD, et al. *In vivo* imaging of light-emitting probes. Journal of Biomedical Optics. 2001;6:432–440. [PubMed]

[23] Schweiger M, Arridge SR, Hiraoka M, Delpy DT. The finite element method for the propagation of light in scattering media: boundary and source conditions. Medical Physics. 1995;22:1779–1792. [PubMed]

[24] Taufiquar AT, Yoon J-R. On uniqueness in refractive index optical tomography. Inverse problems. 2006;22:L1–L5.

[25] Tikhonov AN. Regularization of incorrectly posed problems. Sov. Doklady. 1963;4:1624–1627.

[26] Wang G, Hoffman EA, et al. Development of the first bioluminescent CT scanner. Radiology. 2003;229(P):566.

[27] Wang G, Jaszczak R, Basilion J. Towards molecular imaging. IEEE Trans. Med. Imaging. 2005;24:829–831.

[28] Wang G, Li Y, Jiang M. Uniqueness theorems in bioluminescence tomography. Med. Phys. 2004;31:2289–2299. [PubMed]

[29] Weissleder R, Mahmood U. Molecular imaging. Radiology. 2001;219:316–333. [PubMed]

[30] Weissleder R, Ntziachristos V. Shedding light onto live molecular targets. Nature Medicine. 2003;9:123–128. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |