Computational Engine for a Virtual Tissue Simulator

We have developed a computational platform that simulates light transport in tissue in support of biomedical optics research. Although in its initial stage of development, this platform is being used to answer important questions regarding the detection of tissue changes, and the optimal design and positioning of optical probes to ‘interrogate’ the tissue best. We provide answers to such questions by applying perturbation and midway surface Monte Carlo techniques. Derivation of these methods makes rigorous use of the radiative transport equation which is essential if the methods are to provide accurate solutions for highly complex media such as biological tissue.


Introduction
Computational tools for modeling radiative transport in biological tissues have played a vital role in the development of optical techniques for the diagnosis and therapeutic treatment of tissues.These tools aid in the design of optical probes to detect noninvasively tissue transformations attributed to cancer and other abnormalities.To date, models of tissue have been confined to very simple geometries such as homogeneous and layered media.Recently, however, there is evidence that optical signals provided by multiply scattered light are sensitive to changes in tissue structure and composition on the mesoscopic (0.1-1 mm) spatial scales [FOV + 03, KWR + 03, MYLK05].This realization has driven the need to (a) model tissue with greater spatial refinement, (b) understand the detectability of specific tissue changes, and (c) determine the tissue regions from which the detected light is remitted; i.e., the spatial and angular distribution of the light propagating from source to detector.
We are addressing these needs by developing a virtual tissue simulator (VTS).This computational platform allows the user to specify a probe configuration and define a voxelized tissue representation.A variety of probe configurations will eventually be incorporated.Here, we focus on one consisting of a fiber-optic source and detector at a fixed separation.The voxelized tissue representation can be provided from images provided by histology, CT or MRI.For a specified probe configuration and tissue definition, the user can run a conventional Monte Carlo simulation to model light transport through this system.The VTS also provides perturbation Monte Carlo capabilities that can be used to determine the change in the detected signal due to small changes in tissue structure and/or composition.Finally the VTS incorporates midwaysurface Monte Carlo methods that couple forward and adjoint simulations at an intermediate surface to provide a spatial map of photon propagation from source to detector.Both perturbation and midway-surface Monte Carlo methods provide gains in efficiency and accuracy over conventional Monte Carlo methods.
The long-range goals for the VTS include the ability to solve inverse problems of two types.First, if a particular change in the detected optical signal is detected by a specific probe, we would like to determine the changes in the optical properties of the tissue that produced the optical signal change.We have already demonstrated the use of perturbation and differential Monte Carlo methods [Hay02,HS04,SYHV07] to solve this type of inverse problem and plan to extend this method to systems in which the tissue region is voxelized.Second, if we wish to target a specific tissue region, we would like to determine the probe design characteristics and probe placement that would optimally interrogate the targeted tissue region.Probe parameters such as the radius, numerical aperture or orientation of the source and detector fibers, source-detector (s-d) separation, and relationship of probe to target region can be varied to enhance detection of target tissue changes.Robust forward problem simulations are an essential part of accurate inverse problem solutions.With these ultimate goals in mind, we concentrate in this paper on the forward models for each of these inverse problems.

Monte Carlo Methods for Radiative Transport
The development of a Monte Carlo simulation of light transport in tissue follows from a probability model derived directly from the analytic radiative transport equation.This equation describes the physics of the problem and the probability model defines the probability space needed to solve the problem using Monte Carlo simulations.In the next sections, we present the linkages between the analytic and probabilistic formulations needed to establish the equivalence of these two formulations.The theoretical foundations presented in §2.1 support the Monte Carlo simulations that form the computational engine of the VTS.

Conventional Monte Carlo
The analytic model describing light transport through tissue is the radiative transport equation (RTE).In a closed, bounded subset D of R 3 , the integrodifferential form of the RTE is (1) where Φ(r, Ω) is the photon radiance, with r and Ω representing position and unit direction vectors, respectively.µ t (r) = µ s (r)+µ a (r) is the total attenuation coefficient, µ s (r) is the scattering coefficient, µ a (r) is the absorption coefficient, p(Ω → Ω) is the scattering phase function, and Q(r, Ω) is the internal source function.A unique solution Φ(r, Ω) is assured for all r ∈ D, Ω ∈ S 2 by specifying appropriate conditions on the boundary ∂D.
Detectors are often placed within the tissue system to measure a system response.For example, this system response could consist of reflected and/or transmitted light using of one or more detectors.The quantity describing the system response I in the context of the analytical model is where h is a known 'detector' function, Φ is the solution to Eq. (1) and Γ is the phase space.Suppose, for example, that the amount of energy absorbed within a detector occupying a subregion V of phase space is of interest.In this case, the detector function is defined by h(r, Ω) = µa(r) µt(r) χ V (r, Ω) where and the integral in Eq. (2) measures the amount of photon absorption in the volume V .With this definition of h, I can be estimated as the ratio of the total number of photons absorbed in V and the total launched from the source.
To represent a detector on the surface of the tissue measuring reflectance/ transmittance, the same formalism contained in Eqs.
(2) and (3) can be used.In this case, the volume V is treated as an infinite absorber outside of the tissue whose intersection with Γ is the surface of the detector.The Monte Carlo solution of Eq. ( 1) is captured by a random variable ξ defined on the sample space of all random walks Ω whose expectation is where µ is the analog measure that captures the physical model faithfully and I is given by Eq. (2).Details regarding such constructions may be found in [SG69].

Perturbation Monte Carlo
Once a tissue system of interest is defined within VTS and a Monte Carlo simulation is executed to estimate a desired system response I, the impact on I of changes in the tissue structure and/or composition can quickly be determined using perturbation Monte Carlo (pMC).We briefly review the pMC method next; details can be found in [Hay02, HS04, SYHV07, HSB + 01].
Perturbation Monte Carlo provides radiative transport solutions for multiple systems that can be expressed as a perturbation of a baseline system using a single set of random walks.A set of photon biographies is generated within a baseline tissue system of specified tissue properties.Using pMC, the change to the system response Î = I + ∆I due to perturbations in the optical properties of the baseline system can be determined.This is done by appropriately modifying the random variable used to estimate the reflected light in the baseline Monte Carlo simulation.The pMC method provides estimates of the responses in a 'perturbed' system with a computational cost that is orders of magnitude smaller than that required to run independent Monte Carlo simulations.In addition, the positive correlation between the baseline and perturbed system responses enables pMC to capture small changes ∆I in the system response with a much higher precision than would be obtained with independent simulations.The use of pMC can enable VTS users to determine rapidly the degree to which a diagnostic or therapeutic measurement will be sensitive to changes in tissue properties.
The pMC method can be described within a probability model in terms of the pair of random variables ξ, ξ where ξ is the system response of the baseline tissue system and ξ is the response in the perturbed tissue system.The derivations necessary to make this formulation correct are based on the identity and µ is the analog probability measure based on the baseline optical properties.The measure μ incorporates the analog measure except in the perturbed region, where it uses the optical properties assumed for that region.The Radon-Nikodym derivative (dμ/dµ) expresses how the analog random variable ξ must be modified to produce an unbiased estimator ξ of the optical response in the perturbed system.Explicit formulas derived from Eq. (6) may be found in [Hay02, HSB + 01].
We apply the pMC method to study dysplasia within epithelial tissue.In the initial stages of epithelial tissue dysplasia, cells within the epithelium adjacent to the basal lamina exhibit a larger nucleus to cytoplasm ratio (see Fig. 1) in which the scattering of light is increased by a factor of three [CAM + 03, CFMRK05].Our interest is in determining how the placement of the probe on the tissue surface relative to these dysplastic regions effects our ability to detect these changes.The baseline tissue is modelled as an epithelial region atop a stromal region with an undulating basal lamina interface.The entire tissue is divided into uniform voxels measuring 100 µm × 100 µm in the x-z plane.The voxelized approximation of the undulating interface is shown with white line segments.The optical properties of the two tissue regions are typical of normal epithelial/stromal tissue [CAD + 04]: for the epithelial region, µ a = 0.12/mm, µ s = 2.8/mm with g, the average cosine of the scattering phase function, set to 0.97 and n, the refractive index equal to 1.4.The stromal region optical properties are: µ a = 0.09/mm, µ s = 17.5/mm, g = 0.8 and n = 1.4.The probe configuration consists of a source and detector each with 200µm radius and 0.37 numerical aperture positioned 1 mm apart on the tissue surface.To model the initial stages of dysplastic transformation, one voxel within the epithelial region and positioned atop the two-region interface at x = 0.45 mm (outlined in black) is considered to be 'dysplastic' and assigned a scattering coefficient that exceeds that of the baseline tissue by a factor of 3. The probe is positioned so that the voxel with respect to the x-axis is approximately midway between source and detector: the source is centered at x = 0 mm and the detector at x = 1 mm.This placement is chosen because conventional wisdom suggests that the probe source and detector should straddle the target region.
Fig. 2 displays the detected reflected signal as a function of time, R(t), for both the baseline tissue system and for this system with the dysplastic voxel at x = 0.45 mm atop the two region interface (as shown in Fig. 1).The two plots are visually indistinguishable.This leads us to the conclusion that this particular placement of the probe relative to the voxel exhibiting dysplasia is insensitive to this pre-cancerous tissue transformation.The absence of a significant change in the detected optical signal due to the appearance of the dysplastic voxel can be better understood if one has knowledge of how the light propagates from the source, through the tissue, to the detector.To date, available techniques to provide such information have been based mainly on the diffusion approximation to the radiative transport equation [BOCY97, FZC95, PAEWO95, SHL93].However, the validity of diffusion-based models is compromised when: (a) s-d separations are small or (b) the tissue absorption is comparable to or greater than scattering.In addition, these models cannot, of course, provide transport-theoretic quantitative information.
We have recently developed a novel Monte Carlo method that couples forward and adjoint simulations to generate the spatial distribution of the migration of light from source to detector [HSVed].This map provides valuable information regarding the relationship between a specific probe configuration and placement and the resulting tissue region that can be interrogated.These interrogation density maps are faithful to the radiative transport equation and therefore provide quantitative measures of the contribution of different tissue regions to the optical signal.

Midway Surface Monte Carlo
The magnitude and spatial extent to which an optical probe interrogates the tissue system under examination is of great interest when designing a diagnostic technique.For example, in the case of epithelial dysplasia considered above, we would like to determine the probe position on the tissue surface that would best interrogate the dysplastic voxel.A map that depicts the migration of the light from the source to the detector would aid in the probe placement. 1o generate such a map, each voxel within the tissue representation is treated as a 'target' subvolume.Conventional Monte Carlo methods could be used to select those photon trajectories that have migrated from source to target region to detector either in a forward or adjoint simulation.However, when the source and detector are each small relative to the target subvolume, forward or adjoint simulations used alone engender low statistical signal-to-noise ratios.Such a situation is exceedingly common in biomedical optics.
Our approach to bypass this dilemma, is to break this problem into two components each of which can be determined rapidly through a Monte Carlo simulation.First we determine (a) the probability of photon visitation from the photon source to the target subvolume, P (V ) ('target visitation'), and then (b) the probability of photon detection conditioned by target visitation, P (D|V ), ('detection given target visitation').These two probabilities can be combined using Bayes' Theorem to provide the joint transport probability of 'target visitation and detection': We use a conventional Monte Carlo simulation to determine P (V ) for every voxel in the VTS representation of the tissue.This not only provides an estimate of the P (V ) term in Eq. ( 7), but also a spatially-resolved map of the absorbed or scattered light distribution within the tissue.To determine the second factor P (D|V ), we utilize an adjoint simulation to increase efficiency.This is done by modifying a generalized reciprocity principle that converts P (D|V ) to a coupled forward-adjoint computation at the 'midway' surface of the target subvolume.'Midway' forward-adjoint coupling methods [Cra96, SJH98, UHK01, Wil91] have been successfully applied to increase efficiency in estimating detector responses.In essence, a midway surface between the source and detector is defined such that all detected radiation must pass through this surface.The coupling of a forward and adjoint simulation at this intermediate surface determines the detected response more efficiently, particularly in problems that involve deep penetration.The midway method is made rigorous by utilizing generalized reciprocity theory for transport equations.
We first present the analytical model describing generalized reciprocity and then show how it is modified to allow evaluation of a conditional probability to provide the desired interrogation maps.The intermediate derivations are presented to clarify the final application.

Generalized Reciprocity
Generalized reciprocity establishes the equivalence between the execution of a 'forward' Monte Carlo simulation from a source to detector and an adjoint simulation for 'backward-propagating' photons from the detector (adjoint source) to the source (adjoint detector).This result can be appreciated by first considering the equation adjoint to Eq. (1): where Φ * is the adjoint photon radiance and Q * is any adjoint source function.Let V M be an arbitrary closed, bounded subset of D and ∂V M its surface.If we multiply the radiative transport equation [Eq.( 1)] by Φ * , Eq. ( 8) by Φ, subtract the latter product from the former and integrate the difference over all locations and directions within V M , we get where the variables of integration are suppressed but are understood to be (r, Ω) with the spatial vector r ranging over the volume V M .Using Green's theorem to replace the volume integral on the left side of Eq. ( 9) by a surface integral leads to: where n M is the outward-pointing unit vector normal to ∂V M .Eq. ( 10) is often referred to as the global reciprocity theorem [WE77].Note that if V M = D and the boundary conditions at the air-tissue interface cause the integral on the left hand side to vanish, we then arrive at the 'classical' statement of reciprocity: Eq. ( 11) states that one can obtain the same transport estimates by performing either a forward simulation from the source Q to detector Q * or by performing an adjoint simulation from adjoint source Q * to adjoint detector Q.While Eq. ( 10) is valid generally, it becomes particularly useful when V M encloses either the source or the detector region.The surface of V M , ∂V M , can then be identified as a 'midway' surface between source and detector: every photon that is detected from the source must intersect the midway surface.The function ΦΦ * that occurs in Eq. (10) has been called a 'contributon' response function [Cra96, SJH98, UHK01, Wil91, WE77, SJH99, UH01] and can be used to define a function that characterizes transport from source to detector.If V M encloses the source region, and Q * = 0 in V M , the left hand side of Eq. ( 10) is positive, and equals V M ×S 2 QΦ * , which is the adjoint representation of the detector response.If V M encloses the detector region, and Q = 0 in V M , the left hand side of Eq. ( 10) is negative and equals − V M ×S 2 Q * Φ, which is the forward representation of the detector response.

P (V ∩ D) Maps
We extend this generalized reciprocity theory to determine the conditional probability P (D|V ) for an arbitrary target subvolume V enclosing neither the source nor the detector.Since we are interested not in the interrogation of surfaces but rather of subvolumes within a tissue domain, we slightly modify the midway method to facilitate the estimation of the conditional probability P (D|V ).We launch photons at a physical source Q that propagate until they exit the phase-space.Only photon trajectories that have intersected the target subvolume V contribute to the estimate of P (V ).These 'visiting' photons generate an induced source internal to V that produces a surface source Q ∂V on ∂V.This surface source is then paired with the adjoint radiance on ∂V in a bilinear integration that produces an estimate of P (D|V ).The product of the two probabilities P (V ) and P (D|V ) defines the probability that subvolumes within the phase-space are both visited and detected.We use this product to provide the key quantitative information used to assess and compare the characteristics of potential probe designs.
Let Q V denote the source induced in V by photons launched according to the original optical source function Q(r, Ω).This induced source internal to V, Q V , generates a source on ∂V.If we merely replace the source function Q by the source function Q V (r, Ω) and repeat the derivation that led to Eq. (10), we obtain where the radiance Φ is the solution of the RTE [Eq.( 1)] with the induced source function Q V .Recall, Q * (r, Ω) is a detector function; as such Q * (r, Ω) = 0 unless r is on the boundary of the tissue and Ω points outward .Replacing Q * (r, Ω) by Q * (r, −Ω), therefore, defines an adjoint source pointing into the tissue.This in turn generates an adjoint radiance, Φ * (r, −Ω) inside the tissue.This reverses the direction in the arguments of Q * and Φ * in Eq. ( 12), which then reads since Q * = 0 inside V.The estimation of the right hand side of Eq. ( 13) is performed using an adjoint simulation and provides the detected response due to the induced source Q V , or P (D|V ).
The forward simulation of photons exiting an arbitrary target subvolume V is used to determine P (V ) and is matched at ∂V with the adjoint simulation estimate of P (D|V ).The joint probability of visitation and detection P (V ∩ D) [Eq.( 7)] is formed by the product of these two factors.Fig. 3 displays the resulting P (V ∩ D) map using our baseline epithelium/stroma tissue definition.Each voxel (shown by the superimposed grid in Fig 1) is treated as a target subvolume and the joint probability P (V ∩ D) is determined.Notice that the detected light field does not experience significant lateral dispersion within the epithelial region during its propagation from source to detector.This is due presumably to the high scattering asymmetry coefficient g and low scattering coefficient µ s in that region whose combined effect is to produce minimal lateral dispersion.This behavior below the source and detector within the epithelial region suggests that perturbations within this region, in particular the voxels around positions x = 0 and x = 1, would have the greatest effect on the detected signal.On the other hand, the voxels midway between source and detector above the interface are fairly dark indicative of positions of small perturbative effect.This explains why the dysplastic voxel positioned at x = 0.45 mm did not noticably perturb the detected optical signal.
Based on this analysis, we reposition the probe to place the detector directly above the dysplastic voxel: the source is now centered at x = −0.5 mm and the detector at x = 0.5 mm, while the dysplastic voxel remains atop the two region interface at x = 0.45 mm.Fig. 4 displays R(t) for the baseline tissue system with and without the dysplastic voxel.With the new probe position, the change in the detected signal due to the increased scattering in the dysplastic voxel can now be seen.Notice also the strong correlation between the baseline and perturbed plots.The small variations of the perturbed plot from the baseline capture the effect of the change with much more accuracy than if independent Monte Carlo simulations were used for the baseline and perturbed systems.4. Time-resolved reflectance for the baseline tissue system and for the perturbed system with a dysplastic voxel at lateral position x = 0.45 mm atop the two region interface with the probe source and detector centered at x = −0.5 mm and x = 0.5 mm, respectively.

Summary
We have completed the initial development of the computational engine of a virtual tissue simulator that incorporates efficient forward and adjoint Monte Carlo simulations in a voxelized representation of tissue.This new platform provides the biomedical optics researcher with a means to analyze better the size and location of detectable tissue anomalies and to design and position optical probes to capture these changes effectively.The engine is currently comprised of conventional, perturbation, and midway surface Monte Carlo techniques.In contrast to many other available methods, ours are radiative transport equation rigorous and therefore provide more accurate models that are essential for representing complex media such as biological tissue.
A robust forward model representation is a vital component of an inverse problem-solving capability.By adding differential Monte Carlo to the VTS, we plan next to extend the VTS platform to solve inverse problems.The first step in this direction will enable quantitative determination of optical properties in target regions based on measurements from a given probe design.For example, given baseline and perturbed measurements as shown in Fig. 4, we would like to determine what change in the absorption or scattering properties of the dysplastic voxel caused the perturbed measurement.Second, probe design parameters will be determined that best interrogate specific target regions within the tissue.Generating the change to the P (V ∩ D) maps as a function of probe design parameters will enable inverse solutions that identify the best probe configurations to target selected regions within the tissue.

Fig. 1 .Fig. 2 .
Fig. 1.Schematic of the tissue definition with upper epithelial and lower stromal regions separated by an undulating basal lamina interface.The white superimposed grid identifies the 100 µm×100 µm voxelization.The solid black box atop the interface designates a dysplastic voxel in which the scattering is increased three-fold.

Fig. 3 .
Fig. 3. Interrogation density P (V ∩ D) map of the baseline tissue system with a log-scale gray scale bar.The interface between the epithelial and stromal layers is shown using white line segments.A dysplastic voxel atop the interface at x = 0.45 mm with scattering increased threefold is shown in solid black.
Fig.4.Time-resolved reflectance for the baseline tissue system and for the perturbed system with a dysplastic voxel at lateral position x = 0.45 mm atop the two region interface with the probe source and detector centered at x = −0.5 mm and x = 0.5 mm, respectively.