Sensor optimization for fluorescence optical tomography by experimental design methods

With the increasing importance of molecular imaging fluorescence based methods are continuously gaining impact. In fluorescence optical tomography excitation light is injected into the tissue where the fluorophore converts it to radiation of another wavelength. From the emitted light reaching the boundary the 3-D distribution of the fluorophore is reconstructed. This paper aims at finding the optimal spatial distribution of optodes in order to keep their number (hardware costs) low while gaining maximum information from the target object. The implemented algorithm starts with an arbitrary pool of feasible optodes. The optimal subset is searched by minimizing the mutual information between the different measurements. This goal is reached by subsequently removing those sources and detectors which add the least independent information until a stopping criterion is reached. Mutual information is estimated by calculating the inner products between the rows of the sensitivity matrix i.e. the first derivative of the forward mapping with respect to the optical parameters to be reconstructed. We assembled this matrix with a finite element implementation of the diffusion approximation of light propagation in scattering tissues. When starting with an initial pool of 96 optodes regularly spaced on a cylindrical surface and focusing on different target regions within the cylinder, the algorithm always converged towards physically reasonable optimal sets. Optimal source/detector patterns are be presented graphically and numerically.


Introduction
Fluorescence di usion optical tomography (FDOT) is one of the newer imaging techniques with promising ap-plication potential in medicine.FDOT provides the possibility of functional imaging, i.e it not only visualizes anatomical structures but also provides information about physiological states and processes.FDOT utilizes the ability of so-called uorophores to absorb light in a certain wavelength range and to emit photons at a higher wavelength.e excitation light is injected into the sample through a set of sources.A source can either be in contact with the sample's surface (e.g. a waveguide) or it delivers the light in a contactless manner using collimated or divergent light-beams.e excitation light is sca ered and absorbed while spreading in the tissue.At sites where a uorophore is present and active (e.g.inside a tumor), a part of the absorbed light leads to re-emi ion at another wavelength.is secondary light is again sca ered through the tissue and the part which reaches the boundary can be measured by photo detectors.Due to the di use propagation of the photons in tissue [1] light emerging from the uorescent dye widely spreads before it reaches the boundary.is is in contrast to other established imaging techniques like X-rays where the rays travel through the sample of interest in nearly straight lines.e photon di usion has to be considered in a suitable forward model which is the basis for the reconstruction algorithm that seeks to determine the distribution of the uorophore from boundary measurements.
e reconstructed results usually improve when increasing the number of sensors.However, this is only true up to a certain extent as the di use nature of the photon propagation inherently limits the independence of information of di erent sensors and hence the obtainable resolution.Contrary the computation time and the memory needed for reconstruction increases with the amount of data gathered.e goal is to nd a good compromise between image quality and computational e ort.Graves et al [2] have performed some investigations on how the amount of sources and detectors and their distance, respectively, in uences the reconstruction.Later on the method was extended by Lasser et al [3] who applied it to 360 deg projection tomography. is paper presents a di erent approach to optimize the optode conguration of a uorescence tomography system in the sense that it does not just compare di erent optode con gurations but provides an information measure for every single optode o ering greater exibility, therefore.Furthermore, it will be shown how the optimization can be modi ed such that the reconstruction is focused to a given region of interest.

Forward model
One of the most accurate ways to model light propagation is to utilize Boltzmann's transport equation for kinetic gases.e photons can then be treated like independent gas particles leading to the radiation transfer equation.Unfortunately, the photon intensity in the radiation transfer equation is a eld dependent on the spatial coordinates and the direction (i.e. two angles) into which the photons travel.is leads to a discretization with a huge amount of degrees of freedom and requires extensive computing power and memory.erefore, it is common to use an approximation of the transfer equation known as di usion equation [4].Including a spatially variable urophore concentration c the di usion equation reads together with the boundary condition where φ is the photon density eld, κ = (3(µ a,i (x) + c(x) + µ s )) −1 is the di usion coe cient, µ a and µ s are the intrinsic absorption and reduced sca ering coe cient, respectively, is the uorophore's extinction which relates the urophore concentration to photon absorption.ω is the modulation frequency of the light source, ν is the speed of light and q is the source term.R is a factor to incorporate re ections at the boundary whose normal is denoted by n.Equation 1 can be solved e ciently using a nite element discretization.
For uorescence applications, two di usion equations -one describing the propagation of the excitation photons, another one to describe the emission eld -can be coupled.We prefer to write this in an operator (or matrix-like) notation, where A ex and A em describe the propagation of the excitation and emission eld, respectively.
e operator B is an operator converting the photon density from the excitation wavelength to the emission wavelength at sites where a uorophore is present and serves as a source for the emission eld, thus.It de nition reads whereby Q is the quantum yield [5].
Although more elaborate detector models (e.g.[6]) could be used, in this paper a measurement d is de ned as the number of photons leaving the sample at a certain point x D per unit time: If more than one source and one detector are present, we denote the measurement made with the i-th pair of source/detector by d i .

Sensitivity
In order to solve the inverse problem, i.e. the reconstruction of the distribution of the uorophore's concentration c from measurements on the boundary, it is necessary to know the in uence of a change in the concentration distribution on the measurements.In other words the so-called sensitivity, the derivative of the system 3 with respect to c, is needed.Assuming the existance of the inverses of the operators A ex , A em and B as well as their Fréchet derivatives A ex , A em and B one can derive the expression which relates the change of the concentraton δc to a change in the measurement dd i .
In a nite element context, the discretization of the concentration using piecewiseconstant ansatz functions in equation 6 leads to the Jacobian or sensitivity matrix which is denoted by J in this paper.e element J ij describes the e ect of a concentration change in the j-th nite element on the i-th measurement.
In certain applications it is feasible to operate with di erence measurements.A measurement d 0 is made with a base-line concentration c 0 and a second measurement d 1 is performed a er the concentration distribution has changed to c 1 .If the di erence in concentrations is small, the following linearization can be used for reconstruction: where ∆d is a vector of di erence measurements and ∆c the vector of concentrations in the nite elements.is formulation will be used throughout this paper.However, the ∆ is neglected from now on and we understand all measurement and concentrations as di erences to a base-state.

Redundancy reduction
e optimization approach implemented is based on the idea that the di erent measurements should be as independent as possible, i.e. every measurement should result in new information which can be used for the inverse problem.e method was originally developed by Michelini and Lomax and published by Curtis et al [7] One way to quantify the independence of two measurements is to calculate the inner product and the angle, respectively, between the respective rows of the sensitivity matrix J. en the algorithm has to nd that set of measurements which is closest to an orthogonal set.Let M be the set of all measurement indices, i.e. each element of M uniquely de nes one pair of source/detector.Further, let S i ⊂ M denote the indices of those measurements which are made with the i-th source.e cosine of the angle between two measurements, one made with source i and one made with another source, is given by the term which is a real number from the interval [0, 1].If the measurements are orthogonal, the cosine will be zero and the more they depend on each other, the more the cosine will approach 1.Now, consider the average cosine between all measurements made with source i and those measurements made with another source.is will lead to the expression which is again from [0, 1].If the value is close to 1, the information gained using source i can also be gained using other sources, which means that the source i is redundant.erefore, we call ri the redundancy of source i.
e quantity q i := 1 − r i is used as a quality criterion in the optimization algorithm.An analog expression can be derived for the detectors.is is achieved by replacing Si in [9] by Di , the set of indices belonging to the i-th detector.
e optimization algorithm starts with a set of feasible optodes.It then calculates the quality measure for every source and removes the one with the lowest measure from the optode pool.In the next iteration the same is performed for the detector optodes. is is done until a previously de ned amount of sources and detectors is le in the pool, which is considered to be the optimal measurement con guration.

Focusing
Sometimes it is necessary to focus the optode arrangement so as to reach e.g. higher sensitivity or a higher resolution in a speci ed region than in the rest of the sample.
A rather simple approach is to multiply the sensitivity matrix J with a weighting mask f from the right: whereby J F denotes the focused sensitivity matrix.is resultant matrix will then be used in the optimization algorithm described above.In the simplest case f is a binary vector being one in the region of interest and zero everywhere else.Also smooth variations of the mask are possible.Generally we can assume that 0 ≤ f i ≤ 1.

Results
e optode optimization was performed on a cylinder with a height of 90 mm and a radius of 30 mm, which should mimic a small animal.e values of the optical properties are found in table 1. e uorophore concentration was assumed to be spatially constant and rather low. is should model a cylinder without uorescent dye but with a certain amount of auto-uorescence.Equations and estimates of these parameters can be found in references [5,9,10].e absorption due to autouorescence was chosen to be 10 −2 of the uorophore absorption found in the thesis by Joshi [5].
As initial pool of feasible optodes position a regular grid with 48 source and 48 detector nodes was speci ed. e optodes were arranged in a zig-zag like pa ern on six rings with a spacing of 10 mm (see gure 4).e optimization algorithm needs the desired amount of sources and detectors as stopping-criterion.It was chosen to search for an optimal con guration consisting of eight sources and eight detectors.ree di erent focus regions were chosen to demonstrate the optimization.Region A consists of the voxels in the cylinder slice given by 10 mm < z < 20 mm, region B is another slice de ned by −7.5 mm < z < 7.5 mm and region C is the half-cylinder slice −20 mm < z < −10 mm and x > 0 mm.e optimization was performed on a nite element mesh with approximately 30,000 tetrahedral elements.
As can be seen in the plots 2(a)-2(c), the optimization algorithm tends to concentrate the optodes near the focused region.e source optodes are very close to the region of interest and so are the detectors in 2(c).However, the detector optodes when focusing on larger slices as in 2(a) and 2(b) are more spread.
Figures 3(a)-3(c) showes the rank of the optodes, i.e. the iteration in which they were removed from the pool of feasible optodes.One can see clearly that the optodes far away from the focus region, which is drawn in orange, are removed rst while those near or in the region remain longer in the feasible set.
As a quality measure for the the resultant sensitivity matrix, its singular values and its condition number, respectively, can be used.e largest and smalles singular values together with the ratio between them can be found in 2. e full optode con guration has a rather high ratio of 6 10 12 and is rather ill-conditioned, thus.
e focused designs show a SV ratio which is reduced by a factor of at least 10 6 .Table 2: Singular value (SV) analysis for the full sensitivity matrix and the optimized con gurations with focusing to di erent regions.e table lists the largest and smalles singular value as well as the ratio between them which is a measure of stability for matrix inversion.
e ratio of the design focusing on region A (the uppermost cylinder slice) is higher by a factor of approximately 100 than the other two designs.is is probably due to the rather large spread of the optimal detectors ( gure 2(a)).

Discussion
e location of the sensors and detectors is a critical design parameter for FDOT hardware.A con guration determines the sensitivity of the measurements in a given region and also the obtainable resolution.e method presented herein is based on a simulation model and can be used prior to building hardware, thus.Comparisons of di erent hardware con gurations for uorescence tomography that we are aware of were previously reported by Graves et al [2] and Lasser et al [3].eir approaches were based on a singular value (SV) analysis of the so-called weight matrix.is matrix is essentially the sensitivity matrix used in this paper but was obtained using the normalized Born approximation [11].
In contrast to the previous methods, the optimization algorithm used herein operates on the complete deriva-tive of the di usion approximation given by the system 3. erefore, also the change of the excitation eld due to a perturbation in the uorophore distribution is considered. is is neglected, if a rst-order Born approximation is used.e SV analysis implemented by Graves and Lasser requires a singular value decomposition (SVD) of the sensitivity matrix.is demands a large amount of computation time.As an example, the calculation of a single SVD for a matrix of size 48 × 48 × 30000, takes about 50 minutes on an Intel Core2Duo with two processors.erefore, the SV optimization is limited to 2D applications or to rather simple 3D geometries, which can be modelled with fewer nite elements.e redundancy reduction algorithm does not su er from these limitations as it only needs the sensitivity matrix on which it performs basically inner products.e assemby of the sensitivity matrix as well as the inner products can be parallelized with moderate e ort (the basic principle can be found in reference 12, for example) which allows a further speed-up.
Another advantage of redundancy minimization is provides a quality measure for any single optode rather than comparing complete con gurations.is o ers the possibility to choose a superior con guration rst, which could even be obtained with a completely di erent method, and to optimize the arrangement further through removal of optodes exhibiting a poor quality measure.e optimization algorithm presented works top-down i.e. it starts with a full optode arrangement and iteratively discards the optodes having a low quality measure.Such a strategy is easy to implement but there is no guarantee that it reaches the global minimum. is is due to the fact that the algorithm cannot determine the nal result if in a certain iteration the optode with the worst quality was kept and another one was thrown out instead.To ensure nding the global minimum, all possible optode con gurations needed to be tried which is in the order of 10 18 for the initial con guration presented in here and computationally not feasible, thus.e large spread of the detector optodes as it is visible especially in gure 2(a) may appear counter-intuitive at rst sight.However, still unpublished observations show that the intrinsic absorption µ a,i in uences the optimalcon guration.If this parameter is increased, the optimal optodes cluster more around the focused region.
A drawback of this algorithm is that the desired number of optimal sources and detectors has to be speci ed beforehand.However, other stopping criteria as the minimum resolution or the minimum contrast-to-noise ratio could be preferable.A more dynamic approach would be to monitor a certain quality measure (resolution or CNR) during the removal of optodes from the feasible set.en the algorithm could stop when the quality measure shows a steep decrease, i.e. it would remove only those sources and detectors which do not contribute too much to the quality measure chosen.Further investigations have to be done to include those quality measures into our method.

Conclusion
A quite exible algorithm for optimizing the illumination and detection pa ern in uorescence di usion optical tomography has been presented.e algorithm is based on the maximization of independence between the individual measurements.It has been tested with synthetic data and provided reasonable results.

Figure 1 :
Figure 1: Geometry of the simulation model used for optode optimization together with the initial pool of feasible optodes which are arranged in a zig-zag pa ern on six rings with 10mm spacing.All measures are in mm.

Figure 2 :
Figure 2: e result of the optode optimization for three di erent focus regions which are drawn in orange.

Figure 3 :
Figure 3: e rank of the removed optodes encoded in colour for three di erent focus regions (sketched in orange).Optodes drawn in black or red were eliminated early during the optimization while those in yellow and white remained longer in the pool.e optimal optodes have been omi ed.