A deterministic approach to the adapted optode placement for illumination of highly scattering tissue.

A novel approach is presented for computing optode placements that are adapted to specific geometries and tissue characteristics, e.g., in optical tomography and photodynamic cancer therapy. The method is based on optimal control techniques together with a sparsity-promoting penalty that favors pointwise solutions, yielding both locations and magnitudes of light sources. In contrast to current discrete approaches, the need for specifying an initial set of candidate configurations as well as the exponential increase in complexity with the number of optodes are avoided. This is demonstrated with computational examples from photodynamic therapy.

1. Introduction

In biological and medical sciences, diagnostic and therapeutic instrumentation based on visible, near-infrared or near-ultraviolet light (referred to as optical imaging) is of special interest as it is often non-invasive, the cost for the equipment is moderate, and data acquisition is usually fast (compared to, e.g., MRI). However, the high scattering coefficient of biological tissues in the visible spectrum can be a limiting factor for this technique as it causes photons to propagate in a non-deterministic manner. This in turn makes measurements of superficial structures and the selective illumination of deeper regions significantly more complicated. Due to the stochastic nature of the photon paths, the optimal placement of the optodes is not trivial except for simple regular geometries such as cylinders or spheres.

In recent years, there has thus been increased interest in computation-driven optimization of optical hardware operating in strongly scattering tissues. For example, [1] used a singular value analysis (SVA) of the sensitivity matrix (there called “weight matrix”) to compare measurement setups which differed either in the optode spacing or the measurement type (reflectance vs. transmittance setup). In [2], two optode configurations for a hybrid MRI/DOT measurement device were compared, and the number of singular values above a certain threshold was used as a quality criterion. The work in [3] applied the SVA approach to assess different fluorescence tomography setups in two dimensions, while [4] performed comparisons of three-dimensional setups.

A similar problem is faced in photodynamic therapy (PDT) [5], which is used for dermatological and oncological treatments (e.g., esophageal cancers, especially at a stage where surgical intervention is not indicated). It appears attractive to extend PDT to other carcinomas on epi- or endothelial surfaces, and clinical trials are carried out for, e.g., cervix carcinomas. Some other potential candidates, like mesotheliomas of the thoracic cavity (which are typically difficult to treat), represent a special challenge. A major problem is the design of an appropriate light applicator. In the esophageal cavity, which has a simple geometry, the laser light can be readily applied using a cylindrical scattering device. In contrast, the geometry of the intrathoracic cavity is very complex. The application of a standard esophageal applicator is not recommendable because its area of treatment is small and curved and therefore the illumination would be inhomogeneous. Homogeneity of the irradiation is, however, a crucial design criterion, as inhomogeneities can lead to locally ineffective treatment on the one hand and local overdoses on the other hand. This can have serious consequences, including lethal overdoses [6].

In the past, some designs for flexible light diffusers have been proposed, where a typical solution is to use cylindrical or spherical diffusers in a bag filled with a scattering medium. These are applicable after pneumonectomia, if blood accumulations at the surface of the bags are prevented by continuous rinsing [7-10]. Some regions like the sinus diaphragmaticus are difficult to access and have to be illuminated separately. This can be achieved with wedge-shaped illuminators [11]; however, positioning of these illuminators and homogenization of the illumination is difficult. Another approach is to fill the thorax with a biologically non-hazardous scattering medium (e.g., with intralipid), which can also be used for rinsing to avoid blood accumulation. Typically, a spherical diffuser is used for illumination. However, it is difficult to control the dose rate with this approach. This method has also been applied to cases where lung tissue was not resected [8]. In general, both methods try to achieve homogeneous fluence using real-time dosimetry and manual repositioning of the light diffuser. An interesting alternative consists in textile-based diffusers, where special optical fibers are integrated in a textile. These diffusers are very flexible but suffer from inhomogeneous illumination and a low transmission rate [12, 13]. Recently so called “light blankets” with arrays of cylindrical diffusers [14] or a spirally-wound side-glowing fiber [15] embedded in a bag filled with intralipid were presented. They are easy to fabricate but still show inhomogeneities, especially at the corners. Due to the need for a homogeneous fluence rate, it is of great interest to optimize the placement of these fibers to obtain improved illumination using minimal energy.

The purpose of this work is thus to present a general approach to compute adapted optode locations for different geometries, tissue types, and applications. The method is based on considering this task as an optimal control problem for a partial differential equation describing the diffusion of photons in a strongly scattering medium, where the location of optodes are modeled as a continuous “source field”. The crucial step—first proposed in [16]—is to include a penalty term that favors pointwise solutions. In this way, both locations and magnitudes of the light sources to be placed are obtained in a single step. The main advantage of this approach over previously published—discrete—methods is that no initial maximal or minimal configuration needs to be specified (although an allowable region can be enforced), and that a combinatorial problem with exponential complexity is avoided. In addition, the algorithm is not based on stochastic (e.g., Monte Carlo) methods but is fully deterministic, which eases the verification of the outcome significantly. Finally, the approach is flexible and can incorporate a wide variety of objective criteria (e.g., photon flux over a given boundary section) by changing the target functional. The proposed approach is demonstrated in the context of optimizing the illumination pattern in the photodynamic treatment of intrathoracic cancer.

The rest of this work is organized as follows. In section 2.1, the specific mathematical model for photon diffusion is given. Section 2.2 presents our optimal control framework for optode placement, whose numerical solution is described in section 2.3. The setup of the numerical experiments used for demonstrating our approach can be found in section 3, and the results are presented in section 4. A discussion of the proposed method in section 5 concludes the work.

2. Theory

During photodynamic treatment of cancer, a photosensitizer such as Photofrin is injected intravenously. Afterwards, the cancerogeneous site is illuminated with red to near-infrared light from sources in a diffuser which is applied directly on the region of interest, i.e., in the intrathoracic cavity. The absorption of energy by the photo-activable drug leads to the formation of cytotoxic singlet oxygen, which destroys cancer cells selectively. The challenge is to homogenize the light intensity as both under- and overexposure can lead to ineffective treatment [17].

2.1. Mathematical model

We use the diffusion approximation of the radiative transfer equation to model the steady state of light propagation in a scattering medium [18]. This leads to a stationary elliptic partial differential equation for the photon distribution φ ∈ H1(Ω), The geometry of the object is given by the domain Ω ⊂ ℝ, d ∈ {2,3} being the number of spatial dimensions, with boundary Γ whose outward normal vector is denoted by n⃗. The medium is characterized by the absorption coefficient μ, the reduced scattering coefficient μ′, and the diffusion coefficient . The coefficient ρ models the reflection of a part of the photons at the boundary due to a mismatch in the index of refraction. Finally, the source term q models the light emission of the embedded optodes.

For the optimal control approach, we also require the solution p ∈ H1(Ω) of the adjoint equation for given f ∈ L2(Ω). Both equations should be understood in the weak sense.

2.2. Optode placement optimization

Since optodes act as discrete light sources, the source term can be modeled as for q ∈ ℝ+ and x ∈ Ω, 1 ≤ j ≤ N, where δ denotes the Dirac distribution (i.e., ∫ f dδ (x) = f (0) for all continuous functions f). A straightforward approach for optimizing the placement of the optodes (as was done, e.g., in [19]) would identify a set of M ≫ N possible optode locations x1,..., x and chose the best N locations such that a certain performance criterion J(q) is minimized. The corresponding optimal source magnitudes q would then be computed in a second step.

To avoid the combinatorial complexity of this discrete approach, instead of specifying the optode locations beforehand, we optimize the (distributed) source term q directly while adding a penalty term that promotes sparsity of q, i.e., smallness of its support {x ∈ Ω : q(x) ≠ 0}. This has the added advantage that the number N of optodes need not be specified in advance. Since we are looking for point sources, this requires searching for q in the space of regular Borel measures (which includes the Dirac distribution). Following [20], we are thus led to the optimization problem where ℳ(Ω) is the space of regular Borel measures, i.e., the dual of the space C0(Ω) of continuous functions with compact support on Ω, with norm which reduces to for q ∈ L1(Ω). This is related to the well-known fact that L1 norms promote sparsity in optimization. The penalty parameter α controls the sparsity of the solution: The larger α, the smaller the support of q.

Motivated by the application in PDT, we chose as performance criterion the deviation from a constant illumination z in an observation region ω ⊂ Ω such that where φ| denotes the restriction of φ to ω. Due to the linearity of the forward problem, we can take z = 1 Wm−2 without loss of generality. After optimization, the magnitude of the resultant sources can be linearly scaled to achieve the required illumination z. In addition, we restrict the possible light source locations to a control region ω ⊂ Ω, which does not overlap with the observation region ω (i.e., ω̄ ∩ ω̄ = ∅︀), and enforce non-negativity of the source term q (which represents the optodes). This leads to the following optimization problem: It was shown in [21] that this problem has a solution q* ∈ ℳ(ω), which can be approximated by a sequence of functions q ∈ L2(ω) for γ → ∞ satisfying where p is the solution of (2.2) with right hand side f := φ − z and φ is the solution of (2.1) with right hand side q. Equation (2.4) can be solved using a semismooth Newton method which is superlinearly convergent; see [21]. To globalize the Newton method and closely approximate the solution q* of (2.3), we use a continuation scheme in γ where we iteratively solve the problem for an increasing sequence γ, using the previous solution as initial guess.

2.3. Finite element discretization

The discretization needs to account for the fact that the functions q converge to measures as γ increases. We therefore employ the finite element discretization proposed in [22], where the photon density φ and the adjoint variable p are discretized using piecewise linear elements on a given triangulation T, while the source term q is discretized using linear combinations of Dirac distributions centered at the interior nodes x, 1 ≤ i ≤ N(T), of T: In practice, the number of nodes N(T) will be determined by the need to resolve the geometry of the domain and the required accuracy of the solution of the forward model (2.1). Although further refinement of the triangulation increases the number of possible optode locations, the sparsity-promoting property of the minimized functional discourages placing additional optodes. In fact, it was shown in [22] that for a given discretization of the forward model, the computed sources (for γ → ∞) are optimal among all (non-discretized) measures.

Since the linear finite element basis functions form a nodal basis, the right hand side in the weak formulation of (2.1) for a piecewise linear basis function e becomes i.e., the mass matrix is the identity. Introducing the stiffness matrix A corresponding to (2.1) and the observation mass matrix M with entries M = ∫ e dx, we obtain the discrete optimality system Eliminating q using the last equation and applying a semismooth Newton method, cf. [21], we have to solve for (φ+1, p+1) the block system where D is a diagonal matrix with the entries of the vector d, on the diagonal. It can be shown that the semismooth Newton method has converged once d+1 = d holds. After the final p has been computed, the corresponding control can be obtained from (2.4). The complete procedure is given in Algorithm 1.

3. Materials and methods

The optimization algorithm described in section 2.2 is implemented in Python using the open source finite element library FEniCS [23]. The parameters in Algorithm 1 are set to m* = 34 (such that γ* ≈ 1010) and k* = 20. To model a textile-based diffuser, the material parameters in (2.1) are taken as μ = 10−4 mm−1, μ′ = 10−1 mm−1, and ρ = 0.1992. The influence of the parameter α is illustrated by comparing the results for different values of α specified below.

The meshes for the light diffusers containing the optodes are created with the commercial mesh generator Hypermesh™. To demonstrate the behavior of the optimization algorithm for different geometries, we first consider simple two-dimensional spline models which represent the cross-section of an infinitely long pad. This geometry mimics that of an array of parallel cylindrical diffusers embedded in a scattering substrate. Five single-curved models and four double-curved models with increasing curvature κ were created as shown in Fig. 1. The dimensions correspond approximately to a width of 10 mm and height of 120 mm. In all cases, the region ω in which the illumination should be homogenized are the left and right outer lines (indicated in orange in Fig. 1). The region ω where optodes are allowed to be placed is a single line equidistant from both (indicated by a dashed line in Fig. 1). The meshes for the single-curved models of curvature κ = 5, 10, 20, 40, and 60 consist of 61038, 61789, 67160, 80664, and 105322 finite elements, respectively. The double-curved models of curvature κ = 5, 10, 15, and 20 are comprised of 62349, 70735, 82119, and 104220 finite elements, respectively.

Fig. 1

Two-dimensional model geometries (numbers denote curvature κ).

The photodynamic treatment is simulated by embedding the light diffuser model in the intrapleural space of a realistic three-dimensional human thorax model that is constructed from a stack of CT images. The approximate dimensions are: height 100 mm, width 150 mm, thickness 10 mm. The observation region ω is defined as the outer and inner surface of the model, and ω is an interior manifold equidistant from both (see Fig. 2; ω is indicated in purple). The generated mesh consists of 81770 elements.

Fig. 2

Three-dimensional model. The admissible manifold ω for optodes is indicated in purple.

The results are evaluated quantitatively for different values of the sparsity-controlling parameter α. The coefficient of variation c of the resultant photon density φ over the observation region ω and the number N of sources after the optimization procedure serve as quality measures. For the latter, the nodes in the control region ω satisfying q > 10−16 are counted. We compare the results for α ∈ {0.1, 0.01, 0.001} for the two-dimensional models and α ∈ {0.2, 0.4,...,1.8} for the three-dimensional model.

4. Results

The quantitative results for the two-dimensional geometries are given in Table 1 for the single-curved models and in Table 2 for the double-curved models. As can be seen by comparing the number of active nodes N with the total number of nodes for each model, the algorithm indeed produces discrete sources that can be used as optode positions. The obtained coefficients of variation c indicate that a homogeneous illumination of the desired region is possible at least for α < 0.1, demonstrating the feasibility of the proposed approach. The robustness of the algorithm with respect to geometry is illustrated by the fact that the achieved variations do not depend very much on the curvature. It can also be observed how the penalty parameter α determines the tradeoff between the number of active optodes and the homogeneity of the illumination in the region of interest: larger values of α yield fewer optodes but less homogeneous illumination, again independent of curvature.

Table 1

Results for single-curved models. Shown are the number N of active nodes and the coefficient of variation c of the photon density in the observation domain for different curvatures κ and values of α.

κ	5			10			20			40			60
α	0.1	0.01	0.001	0.1	0.01	0.001	0.1	0.01	0.001	0.1	0.01	0.001	0.1	0.01	0.001
N	22	49	59	18	51	83	15	56	66	20	51	62	22	68	147
cv	2.52e-1	1.82e-2	5.40e-3	2.96e-1	2.07e-2	6.17e-3	1.78e-1	1.96e-2	7.95e-3	1.68e-1	2.53e-2	1.09e-2	1.21e-1	2.02e-2	1.57e-2

Table 2

Results for double-curved models. Shown are the number N of active nodes and the coefficient of variation c of the photon density in the observation domain for different curvatures κ and values of α.

κ	5			10			15			20
α	0.1	0.01	0.001	0.1	0.01	0.001	0.1	0.01	0.001	0.1	0.01	0.001
N	12	50	134	19	40	148	26	49	60	33	77	130
cv	1.65e-1	2.48e-2	1.51e-2	2.03e-1	2.88e-2	2.45e-2	2.24e-1	3.27e-2	2.94e-2	4.92e-1	3.47e-2	3.09e-2

The qualitative behavior of the computed sources for each value of α is shown in Fig. 3(a) and Fig. 3(b) for a representative single-curved (κ = 20) and double-curved model (κ = 15), respectively, where the relative strength of the sources is coded by height. (Note when comparing Tables 1 and 2 with Fig. 3 that neighboring active nodes appear as a single peak and thus can be taken as a single optode.) While for the single-curved model and α = 0.1, the distribution of optodes agrees well with the intuitive choice of equally spaced optodes of approximately equal magnitude, the other values indicate that a better illumination can be achieved with stronger sources towards the tips of the model. It should be pointed out that even in the former case, the number of optodes to be distributed is not obvious. For the double-curved models, the results indicate that optodes should be placed preferentially in regions where the curvature changes.

Fig. 3

Optode positions and relative magnitudes (height-coded) for representative single-curved and double-curved models for three different values of alpha (from top to bottom: α = 0.1, α = 0.01, α = 0.001).

For reference, Fig. 4 shows the corresponding photon densities φ (in Wm−2, normalized to unit mean) plotted along part of the observation region (left line in Fig. 1), illustrating how the parameter α and the model geometry influence the homogeneity of the illumination in this region. As expected, photon fluence shows the most pronounced inhomogeneities close to the borders. In the case of the single curved model, a nearly sinusoidal ripple pattern arises in more than 80 % of the target region, while in the double curved model the ripple is superimposed on a step-profile with the steps located approximately at the zero-crossing points of the curvature. With α = 0.1, the peak–peak fluctuations are still around 40 % of the mean value even far away from the borders, which may be considered as unsatisfactory. However, when decreasing alpha to 0.01 or less, the ripple remains within a few percent, which is sufficient, especially when comparing this value to other sources of fluctuations of the irradiation such as local absorption changes by tissue inhomogeneities, bleeding, or inhomogeneities of the distribution of the photosensitizer.

Fig. 4

Photon densities φ (in Wm−2, normalized to unit mean) plotted along part of the observation region (left line in Fig. 1) for representative single-curved and double-curved models for three different values of alpha (from top to bottom: α = 0.1, α = 0.01, α = 0.001).

The quantitative results for the three-dimensional model are shown in Table 3. For α = 1.8, no controls are placed and thus the photon density is zero. This is consistent with the theory, which predicts that there is a threshold value for α above which the optimal control is identically zero; cf. [22, Proposition 2.2].

Table 3

Results for three-dimensional model. Shown are the number N of active nodes and the coefficient of variation c of the photon density in the observation domain for different values of α.

α	1.8	1.6	1.4	1.2	1.0	0.8	0.6	0.4	0.2
N	0	12	150	250	333	409	498	637	884
cv	—	1.85e+0	5.64e-1	3.59e-1	2.65e-1	2.04e-1	1.56e-1	1.13e-1	6.72e-2

Figure 5 shows location and magnitude (color coded) of the computed optodes and the corresponding photon densities (in Wm−2, normalized to unit mean) for α = 1.2, α = 0.8, and α = 0.4. Due to the nonuniform curvature of the model, a homogeneous illumination is harder to achieve than in the two-dimensional case, especially at the borders of the target region. However, for α < 1.2, the inhomogeneities in the interior are usually within 10 %, and the few hot spots of 30 % would still be acceptable. Although of course the specific placement may be difficult to realize in practice, the qualitative distribution can be useful information in the initial design process.

Fig. 5

Optode positions and magnitudes (left) and photon densities (right; in Wm−2, normalized to unit mean) for the three-dimensional model and three different values of α.

5. Discussion

The proposed approach is able to generate reasonable optode configurations adapted to specific geometries, even in situations (such as complex three-dimensional models) where optimal setups are not intuitively obvious. Our method also yields relative strengths of the optodes to be placed, which would otherwise have to be computed in a separate step. Furthermore, the algorithm is deterministic and does not require a-priori knowledge such as an initial set of candidate locations or the number of optodes after optimization, which on the contrary is provided by our approach. The method can be used as a tool during the initial design process to estimate the number of sources required as well as their location and relative strengths.

By formulating the optode placement problem as a continuous optimization problem, the combinatorial complexity inherent in discrete approaches is avoided. This is critical for achieving an efficient optimization technique and—to our knowledge—has not been presented before in the context of diffuse optical imaging. As an example, our Python implementation required about three minutes on a MacBook Pro (2.16 GHz Intel Core2 Duo with 2 GByte RAM) for the single-curved model with κ = 5. Our approach could therefore also be used in an interactive setting, where the engineer will adapt design parameters, such as the optical coefficients of the diffuser, based on the outcome of an optimization run.

While the number of desired optodes is correlated with the penalty parameter α, it is not directly controllable. This drawback is analogous to the problem of finding the “best” regularization parameter in image reconstruction (e.g., for diffuse optical tomography), where typically the determination of the parameter is left to the user or is based on heuristics. Certainly, one could think about finding a good parameter through successive optimization runs, e.g., with decreasing values of α if the user specified an upper bound on the number of optodes.

The achieved results are satisfactory from the mathematical point of view; but of course they should also be discussed in an engineering context. In particular, it may be difficult to place many sources in a more or less irregular pattern. The number of optodes depends on the required uniformity of the surface fluence. A reasonable value in practice would be a CV of 0.05. Table 1 shows CVs (for single-curved pads) below 0.03 for around 50 optodes, but up to 0.25 for less than 25 optodes. The numbers for the double-curved pad are only slightly greater. This means that the between 40 and 50 required sources can be expected. In practice, such a design can be approximated comparatively easily with parallely arranged cylindrical polymer diffusers of sufficiently small radius, which are fed by individual optical fibers. Instead of a fixed grid of diffusers, one can imagine dense bundles of uniformly spaced diffusers where only those close to the optimal positions are connected to the laser source. This would allow a very flexible use and adaptive homogenization of the fluence dependent on the individual anatomical situation (e.g., curvature), which is certainly desirable in the context of a personalized optimization. Such a concept can be realized by using fiberoptic switches with many channels. In three dimensions, the sources may be fiber-coupled spherical diffusers or simply open-ended fibers. Due to the higher number of potential positions (637 for a CV of 0.11, see Tab. 3), the construction of a flexible structure here may be difficult, and pre-fabricated pads that are adapted to a certain anatomical target geometry appear more realistic than a truly adaptive system.

Although a rigorous sensitivity analysis has not yet been carried out, our experience indicates that the computed photon density distributions are relatively robust to small perturbations of the optode locations and magnitudes. Similarly, we did not observe significant changes in the results due to small random perturbations of the optical parameters. This can be attributed to the linearity and the strong diffusivity of the model (2.1). Such robustness is very important for practical implementations because it means that the result is not very sensitive to manufacturing tolerances.

One of the main advantages of the optimal control approach is its flexibility. For example, it is straightforward to extend the underlying model to include, e.g., inhomogeneous material properties or to replace the diffusion approximation by a more complicated model such as the radiative transfer equation. It is also possible to consider different objective criteria such as the photon flux through a given (boundary or internal) surface by changing the functional J(q). In principle, the approach can be applied to the problem of optimal experiment design for optical tomography, if the objective J(q) is based on a suitable sensitivity term. However, this extension of our method is subject to future work.

Algorithm 1 Semismooth Newton method with continuation

1: for m = 1,...,m* do

2:  set γ = 2(m−1), φ0 = p0 = d0 = 0

3:  for k = 0,...,k* do

4:    solve (2.5) for φ k+1, pk+1

5:    compute dk+1 from (2.6)

6:    if dk+1 = dk then

7:      set q(m) = γ min(0, pk+1 |ωq + α)

8:      break

9:    end if

10:  end for

11:end for