Sensitivity Matrix Implementation in the Monte Carlo Package MBioICFO for the Simulation of Photon Migration in Tissues

! ! Master Erasmus Mundus in Photonics Engineering, Nanophotonics and Biophotonics. Europhotonics MASTER THESIS WORK Sensitivity Matrix Implementation in the Monte Carlo Package MBioICFO for the Simulation of Photon Migration in Tissues Niklas Gerdes Supervised by Prof. Turgut Durduran, (ICFO) Co - supervised by Prof. Uli Lemmer, (KIT) Presented on 23rd of October, 2015 Registered at I hereby declare that the present thesis is original work written by me alone and that I have indicated completely and precisely all aids used as well as all citations, whether changed or unchanged, of other theses and publications. Ich versichere hiermit, die vorliegende Arbeit selbstst¨ andig angefertigt, alle benutzten Hilfsmittel vollst¨ andig und genau angegeben und alles kenntlich gemacht zu haben, was aus Arbeiten anderer unver¨ andert oder mit Ab¨ anderungen entnommen wurde. Niklas Gerdes Barcelona, 15.10.2015 Abstract An extension of the Monte Carlo method in diffuse optics was developed. Diffuse optical technology measures light absorption and scattering in human tissue. A sensitivity matrix has to be constructed to obtain structural information. It contains the sensitivities of the detected signal to absorption or scattering changes in different regions of the tissue. The existing MBioICFO simulation package was extended to allow the construction of the sensitivity matrix for perturbations in absorption. As opposed to prior implementations, the sensitivity matrix is determined in one simulation run and for arbitrary geometries. The new method was verified with analytical solutions for homogeneous media with infinite and semiinfinite boundary conditions. The method also enables determination of the sensitivity matrix for different detection times and in geometries obtained from MRI measurements. The program has shown appropriate computational efficiency with acceptable runtimes. iii iv ABSTRACT Acknowledgments First and foremost I want to thank Prof. Turgut Durduran for giving me the opportunity to work on this project in his research group. I owe my progress to his patience with someone completely new to the field of diffuse optics. Both academically and personally, I am grateful for the experience I have had in the medical optics group at ICFO. The research group provided me a supportive and motivational work environment. I would like to point out the supervision and guidance I received from Dr. Johannes Johansson throughout my project. Nicolas Mateos also helped me a lot during the initial setup of the program. Besides, I want to thank the whole ICFO community for making my time here so enjoyable. The numerous extraordinary personalities that one constantly gets to meet at ICFO provide the perfect environment for academic and personal growth. I thank Prof. Uli Lemmer at KIT for his willingness to co-supervise this work. I also want to express my deep gratitude to everyone making the Erasmus Mundus Europhotonics program possible. In the two years I have been part of this program, I have studied and lived in three different countries, built relationships with people from all over the world and expanded my understanding of science. This amazing opportunity I have been given inspires me for all my future work. I hope I can give something back. I could not have done any of this without the people closest to me. You know who you are. You are constantly with me in my thoughts. v vi ACKNOWLEDGMENTS Contents Abstract iii Acknowledgments v Contents 1 1 Introduction 3 2 Theoretical Background 7 2.1 DiffuseOptics............................ 7 2.2 Image Reconstruction and the Jacobian . . . . . . . . . . . . . . 12 3 The Monte Carlo Method 19 3.1 Implementation........................... 19 3.2 Input................................. 24 3.3 Output................................ 25 3.4 Variations in the Monte Carlo Method . . . . . . . . . . . . . . . 27 4 Numerical Solution of Jacobian 29 4.1 Theoretical Approach . . . . . . . . . . . . . . . . . . . . . . . . 29 4.2 Implementation........................... 33 5 Evaluation 37 5.1 The Jacobian in the Infinite Medium . . . . . . . . . . . . . . . . 37 5.2 The Jacobian in the Semi-infinite Medium . . . . . . . . . . . . . 41 5.3 The Jacobian at Different Detection Times . . . . . . . . . . . . . 44 1 2CONTENTS 5.4 The Jacobian in MRI of Human Head . . . . . . . . . . . . . . . 46 6 Conclusion 51 Bibliography 53 Chapter 1 Introduction There are various techniques to image living tissue and ”see” the human body in a way invisible to our bare eyes. In diffuse optics, near-infrared (NIR) light is used to obtain physiological information about biological tissues [8]. Photons in this spectral range (∼650-900 nm) experience low absorption in water and hemoglobin and can therefore travel deep into tissue [27]. Starting from a light source on the surface, photons propagate through the tissue, are scattered multiple times and are ultimately either absorbed, leave the tissue, or are detected by a detector placed at short distance (∼mm-cm) from the source. On their way from source to detector, these photons have experienced absorption, scattering and phase shifts from moving scatterers, all of which influence the detected signal. Hence, analysis of the signal provides information about tissue absorption, scattering and flow of scatterers. Typical absorption lengths after which a photon is absorbed are several centimeters while scattering occurs at distances less than a millimeter. Main absorbers in biological tissue are water, melanin and hemoglobin. Scattering predominantly happens at cell nuclei and mitochondria, since the refractive index difference to the surrounding water or lipid is large. Diffuse optical technology can measure changes in absorption and scattering and therefore changes in the concentrations of absorbers and scatterers. This can be relevant, for example, to measure blood oxygenation [2]. Since oxyand deoxyhemoglobin have different absorption spectra, measuring absorption at several wavelengths allows to determine their con3 10 CHAPTER 2. THEORETICAL BACKGROUND fluence rate is obtained in the following way: Φ(~r, t) = Z4π L(~r, ˆ Ω, t)dΩ(2.4) To summarize, the approximations made to derive the diffusion equation are the following: a much bigger reduced scattering coefficient than absorption coefficient (a factor of at least 10), an isotropic source, slow temporal variations, photon propagation distances larger than the transport mean-free path, scattering angles independent of initial photon direction and observation points far from sources and boundaries. The usual way to solve the diffusion equation is to employ the Green’s function. This is the function that gives the fluence rate for an infinitely short-pulsed point source S(~r, t) = δ(~r, t). The form of the Green’s function depends on the imposed boundary conditions. For an infinite, homogeneous medium, the Green’s function in the time-domain is: G(~r, t) = v (4πDvt)3/2exp −r2 4Dvt −µavt(2.5) This function describes the broadening due to scattering, except for the exp(−µavt) term which represents absorption. To find solutions for any arbitrary source S(~r0, t0), we take the convolution of the Green’s function and the respective source term: Φ(~r, t) = Zt 0Z∞ 0 G(~r, t,~r0, t0)S(~r0, t0)d~r0dt0(2.6) The Green’s function in the continuous wave regime is found by eliminating the time-derivative term on the left-hand side of the diffusion equation 2.2. The Green’s function then reads: G(~r,~r0) = 1 4π|~r −~r0|exp(−k|~r −~r0|)(2.7) In this equation, kis defined as k≡pµa/D. Obviously, one is interested in geometries that model biological tissues. These include boundaries which means that appropriate boundary conditions have to be found. One already very useful geometry is the semi-infinite medium, assuming 2.1. DIFFUSE OPTICS 11 one boundary plane between turbid medium and for example air. To find boundary conditions, one considers that photons leaving the tissue will never re-enter it again. Hence, all incoming radiance is due to Fresnel reflections at the interface. From this consideration, one obtains the partial-flux boundary condition which is exact but cumbersome to handle. A more practical approximation of it is called the extrapolated-zero boundary condition: Φ(z=−zb) = 0 (2.8) It is derived from a Taylor expansion of the fluence rate in the partial-flux boundary condition. It states that the fluence becomes zero at a plane parallel to the interface outside the turbid medium at z=−zb, where the z-direction is perpendicular to the boundary plane with positive z-values inside the turbid medium. The distance zb= 2ltr(1 + Reff )/3(1 −Reff )depends on ltr and Reff . Collimated beam sources at the surface are approximated as isotropic sources at a depth z=ltr, when the photon paths are completely randomized. The coefficient Reff depends on the refractive indices of turbid and adjacent medium. If there is no index mismatch, Reff = 0. For a refractive-index-mismatched boundary with n=nin/nout, the effective reflection coefficient can be approximated as Reff ≈ −1.440n−2+ 0.710n−1+ 0.668 + 0.00636n. From electrostatics, we know the method of image charges to fulfill different boundary conditions. In analogy to that, we can use image sources to fulfill the extrapolated-zero boundary condition. The idea is to introduce a second, negative source outside the turbid medium. Placing a negative point source at position zs=−(2zb+ltr), the superposition of the two point sources makes the fluence vanish in the desired plane. The two point sources are themselves infinite medium solutions, but their superposition is the semi-infinite medium solution. Therefore, the Green’s function in the continuous wave regime for the semi-infinite medium is: G(~r, r1, rb) = 1 4πexp(−kr1) r1 −exp(−krb) rb(2.9) This is obviously merely the sum of two point sources as depicted in equation 2.7. We assume that both sources lie in the origin of the xy-plane. Accordingly, polar 12 CHAPTER 2. THEORETICAL BACKGROUND coordinates are a good choice of coordinate system with ρ2=x2+y2and the source position variables r1and rbare: r1=p(z−ltr)2+ρ2(2.10) rb=p(z+ 2zb+ltr)2+ρ2(2.11) Measurements can be carried out with source-detector pairs either in transmission or reflection geometry. In reflection geometry, light is detected on the same surface where it was injected by the source, some distance ρaway. In transmission geometry, on the other hand, the detector is placed on a surface parallel to the one where the light is injected. For the diffusion equation to be valid, source-detector separation should correspond to at least three times the transport mean-free path ltr. Tissue measurements are done with three different light sources: continuous wave (CW), time pulsed for time-resolved spectroscopy (TRS) and intensity modulated to measure in the frequency-domain (FD). CW sources provide a constant intensity and information is obtained from measuring the drop in intensity some distance ρaway from the source. They are simple and easy to handle, but µaand µ0 scannot be determined simultaneously. Time pulsed sources emit a short light pulse of the order of less than 100 ps. Propagating through the medium, the pulse broadens and detected photons provide information about different locations in the medium depending on when they are detected. In TRS, both µaand µ0 scan be determined simultaneously. Intensity modulated sources, finally, contain the same information content as time pulsed sources. Thus, µaand µ0 scan also be found in a single measurement. The light intensity is modulated sinusoidally at angular frequencies of between 100 MHz to 1GHz. One determines the optical properties by recording the amplitude and phase changes. Time pulsed light and intensity modulated light are simply related via a Fourier transform. 2.2 Image Reconstruction and the Jacobian The general purpose in tomography is to reconstruct a three-dimensional image of the object under investigation. In diffuse optics, there is a mathematical formalism 2.2. IMAGE RECONSTRUCTION AND THE JACOBIAN 13 to reconstruct an image of the optical properties from fluence rate measurements of various source-detector pairs [2]. When traveling through tissue, photons move on a random walk with mean step size ltr. The step size as well as the scattering angles obey probability distribution functions. So their actual values at a specific scattering event cannot be predicted exactly. A visualization of how a photon moves through a turbid medium can be seen in figure 2.2. Source 1/ ' s 1/ a Absorber Scatterer Figure 2.2: Illustration of photon random walk for a detected and an absorbed photon. Note that typically 1/µa1/µ0 sfor turbid media in diffuse optics. In figure 2.3, we can see the regions in the turbid medium that most detected photons have passed on their way from source to detector. For a source-detector pair in reflection geometry and with a certain distance ρ, this region has a banana-like shape. For zero source-detector separation, it resembles a drop. Naturally, the optical properties of the tissue in these regions will affect the detected signal more than those in parts of the tissue that few or no detected photons have visited. To obtain structural information about the tissue from the measured signal, one needs to know how perturbations in the optical properties at different locations influence the fluence rate at the detector. In the following, I will present the theoretical framework of frequency-domain image reconstruction in an infinite 14 CHAPTER 2. THEORETICAL BACKGROUND medium. Source Detector ρ Figure 2.3: Illustration of path distributions of detected photons We assume that the tissue we want to image has perturbations of the absorption coefficient µaonly (equation 2.12). For contrast in the reduced scattering coefficient µ0 s, the derivation is done in a similar fashion. µa(~r) = µa0+δµa(~r)(2.12) This means that there is a spatially dependent perturbation δµa(~r)additional to the constant background absorption coefficient µa0. We assume the perturbation to be small compared to the background, δµa(~r)µa0. The measured fluence rate U(~r)contains a part U0(~r)caused by the background and a part Usc(~r)caused by the perturbation. The Born approach formulates the fluence rate as U(~r) = U0(~r) + Usc(~r), whereas the Rytov approach states it as U(~r) = U0(~r) exp(Usc(~r)) [12]. We will focus on the Born approach. In the forward problem, the fluence rate change due to the perturbation in absorption is calculated. For image reconstruction from fluence rate measurements, however, the inverse problem has to be solved. That means that from the change of the fluence rate, the perturbation δµa(~r)is determined. We approximate the fluence rate as a 2.2. IMAGE RECONSTRUCTION AND THE JACOBIAN 15 first-order Taylor expansion: U(~r) = U0(~r) + ∂U0(~r) ∂µa δµa(2.13) For practical reasons and since spatial resolution is low, the object to be imaged is divided into volume elements, so-called voxels. The fluence rate at the detector has a different sensitivity ∂U0/∂µafor each voxel. The matrix [W] containing the sensitivities of all voxels is called sensitivity matrix or Jacobian. It is the link between the perturbations in absorption and the fluence rate at the detector: [Usc(~r)] = [W][δµa(~r)] (2.14) Equation 2.14 written out explicitly results in: Usc(~rsi, ~rdi) = NV X j Wijδµa(~rj)(2.15) The index i refers to the source-detector pair while the index j refers to the respective voxel inside the medium. The position vectors ~rsand ~rddenote the location of source and detector, respectively, NVis the number of voxels. Determining the Jacobian numerically for any kind of asymmetric geometry using Monte Carlo simulations is the goal of this work. To verify the results though, we can set up an analytical expression of the Jacobian for a homogeneous medium. The approach is to plug the expressions for the fluence rate according to Born U(~r) = U0(~r) + Usc(~r)and expression 2.12 into the diffusion equation 2.2. This leads to: (∇2−k2)Usc(~r) = δµa(~r) DU(~r)(2.16) Following the usual Green’s function approach, we can solve equation 2.16 by taking the convolution of the right-hand side and the corresponding Green’s function: Usc(~rs,~rd) = Z−δµa(~r) DG(~rd,~r)U(~r,~rs)d3r(2.17) Once again discretizing the object into voxels, we can write the integral as a sum 16 CHAPTER 2. THEORETICAL BACKGROUND just as shown in equation 2.15. The sensitivities Wij are then: Wij =∂U0 ∂µaij =−∆V DG(~rdi,~rj)U0(~rj,~rsi)(2.18) The voxel size is denoted by ∆V. The number of rows of matrix [W] equals the number of source-detector pairs, while the number of columns is the number of voxels. For a point-source, it should be noted that in equation 2.18, the expressions G(~rdi,~rj)and U0(~rj,~rsi)are represented by the same Green’s function, with the only difference that they depend on different position coordinates. This solution is illustrated in figure 2.4. First, we take the fluence rate of the source at voxel position ~rj,U0(~rj,~rsi). In that voxel position, we assume the source −δµa(~r) DU0(~rj,~rsi) whose fluence rate is described by G(~rdi,~rj). Integration over space leads us to the final expression in equation 2.17. Figure 2.4: Illustration of analytical expression of Jacobian. Eventually, however, one is interested in determining the image of absorption contrast of the tissue, represented by the perturbation vector [δµa(~r)], from the fluence rate measurements [Usc(~r)]. To do that, we need to invert [W]: [δµa(~r)] = [W]−1[Usc(~r)] (2.19) Inverting [W] is non-trivial and computationally expensive. Among other techniques, singular-value decomposition is often used for this purpose. In this work, 2.2. IMAGE RECONSTRUCTION AND THE JACOBIAN 17 the focus is on numerical determination of the Jacobian [W], not on its inversion. Constructing the proper Jacobian is one of the main steps in the tomography problem. 18 CHAPTER 2. THEORETICAL BACKGROUND Chapter 3 The Monte Carlo Method In diffuse optics, one usually employs the diffusion equation 2.2 as presented in section 2.1 to solve any problem analytically. But it was derived from the radiation transport equation 2.1 based on several assumptions, mainly that µaµ0 s and that the fluence rate is observed far away from sources and boundaries. We can, however, also solve the radiation transport equation numerically using the Monte Carlo method. Monte Carlo simulations are accurate and practically only limited by computational speed. In the following, I will give a thorough presentation of the Monte Carlo method as nowadays used in biomedical optics. I will focus on the open source code MBioICFO [21] that I used and modified during this work. It is based on the widely used Monte Carlo for Multi-Layered media (MCML) approach as presented by Wang et al [26]. MBioICFO is implemented in the object-oriented programming language C++. Since many different Monte Carlo approaches to simulate light propagation in tissue have been published over the years, I will point out differences between them and especially different ways of optimizing computational efficiency. 3.1 Implementation Photon transport in highly scattering media is governed by random processes. Such are the path length before the photon is scattered or absorbed, the direction 19 26 CHAPTER 3. THE MONTE CARLO METHOD allows the calculation of physical quantities for different optical properties of the material without running the simulation more than once. As an example, to calculate the CW intensity at the detector for different absorption coefficients of the materials, we run one simulation with an absorption coefficient of zero of all the materials. With the information of the pathlengths through the materials provided by the history file, we can calculate the weights of the individual photon packets detected with the Beer-Lambert’s law: I= Np X i=1 Wi= Np X i=1 exp − Nm X j=1 µa,jsi,j(3.8) So in the postprocessing, we can directly substitute different absorption coefficients µa,j for any of the Nmmaterials without losing time waiting for more simulations to finish. The pathlengths si,j are all given in the history file. Moreover, we can use the history file to compute the electric field autocorrelation function G1(τ)in the detector: G1(τ) = Np X i=1 exp − Nm X j=1 µa,jsi,j·exp − Nm X j=1 α 3µ0 s,jsi,jh∆r2 j(τ)i2πnj λ02 (3.9) This includes the factor αaccounting for the ratio of scattering that happens at moving scatterers to scattering at static ones. One often assumes it to be unity. The index of refraction njof the respective material as well as the vacuum laser wavelength λ0are included. Again, one single simulation run is sufficient to compute the electric field autocorrelation function for different absorption coefficients µa,j and different mean square displacements h∆r2 j(τ)iof the scatterers. Assuming Brownian motion for example with h∆r2 j(τ)i= 6Dbτ, different diffusion coefficients Dbcan be substituted into equation 3.9 for all the Nmdifferent materials. Alternatively, MBioICFO provides the option to generate an autocorrelation file for the entire geometry, so that one autocorrelation function with various delay times is determined for every voxel. This, however, leads to huge output files which might be difficult to handle. Usually, especially if only the autocorrelation 3.4. VARIATIONS IN THE MONTE CARLO METHOD 27 in the detector is of interest, one should use the history file. Still, the autocorrelation file can be useful to determine the sensitivity S(τ,~r0)at voxel position ~r0and source and detector at positions ~rsand ~rd, respectively: S(τ,~r0) = G1(τ,~rs,~r0)G1(τ,~rd,~r0)(3.10) Unfortunately, two simulations have to be run and hence two large files have to be handled. This is unpractical and ultimately we want to be able to do this in one simulation run. 3.4 Variations in the Monte Carlo Method Most Monte Carlo simulation packages designed for biomedical optics are based on the approach outlined above. Most efforts on expanding and improving the algorithms concentrate on making the method computationally more efficient. One very effective approach is to use information from only one simulation and modify optical properties in the postprocessing. This is what I explained in the previous section with the history file. Zaccanti et al [19] have demonstrated a similar approach. They report a method in which the locations of all scattering events are recorded during the simulation. The temporal response when scattering or absorbing perturbations are introduced is then evaluated in the postprocessing using two scaling relationships, one both for absorption and scattering. So again, information from only one simulation run is needed. Boas et al [3] have presented a simulation package ”tMCimg”, in which spatially varying optical properties in 3D media can be introduced to solve the forward problem. Fang et al [9] have extended this code to a package called ”Monte Carlo eXtreme” (MCX) for parallel computing to achieve shorter runtimes. Other groups also report an acceleration of computation times by a factor of up to 102103by using graphics processing unit (GPU) based Monte Carlo implementations [1,18]. An alternative to the voxelized model is a mesh-based geometry as shown by Margallo-Balb´ as et al [16] in their code ”TriMC3D”. To improve computational efficiency, they use a geometry based on a set of triangle meshes structured with 28 CHAPTER 3. THE MONTE CARLO METHOD a space partitioning scheme. Wang et al suggested using the Monte Carlo method in conjunction with diffusion theory in a hybrid model [25]. This method combines the advantage of the Monte Carlo method, accuracy, and of diffusion theory, computational efficiency, to reach up to 100 times faster computation times than the conventional Monte Carlo approach. In the hybrid model, the Monte Carlo method is only used in regions where the diffusion approximation does not hold. Besides, there exist various techniques in the simulation of photon propagation and detection to reach better computational efficiency. One way of increasing efficiency, for example, is to split photons during their random walk [24]. This makes sense for large source-detector separations so that more split photons actually reach the detector. To conserve energy, the weight is divided on the split photons. A very similar technique is called forced detection as reported by Churmakov et al [5]. This method calculates the small probability that a photon packet goes directly from a scattering event to the detector and a weight is detected proportional to that probability. Various simulation packages also make use of symmetry in the simulated geometry to reduce computation time. Axial symmetry around the light source can be used for example [15]. Some groups have also addressed computing the Jacobian within the Monte Carlo method. The standard approach is to run two simulations to solve the adjoint problem, with the source at the detector position for the second simulation run. The fluence rates of both simulations in each voxel are then multiplied, giving the elements of the Jacobian. Alternatively, axial symmetry can be exploited by introducing absorption perturbations successively in voxels on a radial line and recording the thereby created change in fluence rate. Zaccanti et al [19] developed a code that records the location of all scattering events to plot the scattering density. While this comes close to the Jacobian, it does not properly represent it. I am aiming to implement a numerical solution that computes the correct Jacobian of any asymmetric geometry in one simulation run. One could then use this Jacobian in the inverse problem. Chapter 4 Numerical Solution of Jacobian For finding an implementation of the Monte Carlo method that numerically solves the Jacobian, we face two challenges. First, we have to theoretically develop an approach of how to correctly calculate the Jacobian within the existing Monte Carlo method. How can we use the information about photon transport in tissue provided by the simulation? Next, we have to ensure the developed approach is computationally efficient. This is the major limiting factor of Monte Carlo simulations. Both problems have to be solved to end up with a useful code. In the following, I devote one section to each of the two problems. 4.1 Theoretical Approach As explained in section 2.2, each element of the Jacobian gives the differential of the fluence rate at the detector with respect to a perturbation in absorption in a specific voxel of the geometry [8]. So it is a measure of how much the fluence rate at the detector changes when the absorption coefficient in one or several voxels changes. Naturally, the Jacobian elements for voxels that many photons pass through on their way from source to detector will have a higher absolute value than those few photons pass through. The Monte Carlo method simulates the path of all photons through the tissue. This photon path information is valuable for constructing the Jacobian. For a photon path length lithrough voxel iof absorption coefficient µa,i, the prob29 30 CHAPTER 4. NUMERICAL SOLUTION OF JACOBIAN ability of absorption according to Beer-Lambert’s law is P(absorption in voxel i)=1−exp(−µa,i ·li)(4.1) The probability of the photon to be absorbed at any point on the path from source to detector is P(absorption anywhere)=1−exp(− NV X i=1 µa,i ·li)(4.2) The path length liwill be zero for most of all NVvoxels since a photon packet usually only visits a fraction of all voxels in the geometry. The sensitivity of the fluence rate at the detector to absorption in voxel iis then for one photon: Ji=1−exp(−µa,i ·li) 1−exp(−PNV i=1 µa,i ·li)(4.3) The task is to find a way of implementing the computation of this sensitivity into the Monte Carlo method which uses the concept of weights to account for absorption and Beer-Lambert’s law. Let us start by considering the scattering density ns(~ri, t)as already introduced by Zaccanti et al [19]. It gives the density of scattering events detected photons experienced in voxel iat position ~riwhen detected at time t. It is normalized by the total number of photons launched at the source: ns(~ri, t) = number of scattering events in voxel i total number of photons launched (4.4) With the mean step size of 1/µtbetween two scattering events, we can approximate the average path length liin each voxel as: li=ns(~ri, t) µt (4.5) Again using Beer-Lambert’s law, we can write the relative decay of the number of photons Iiin voxel ias: Ii=Ii−1exp(−µa·ns(~ri, t) µt )(4.6) 4.1. THEORETICAL APPROACH 31 This leads us to formulate the change of the number of photons when an absorption perturbation is introduced: ∆Ii=Ii−1(exp(−(µa+ ∆µa)ns(~ri, t) µt )−exp(−µa·ns(~ri, t) µt )) (4.7) Here, we assume that µaµsand ∆µa< µaso that the change in the total attenuation coefficient µt=µa+µscan be neglected. For small perturbations ∆µawe can make use of the first order Maclaurin series ex≈1 + x: ∆Ii=Ii−1(µa·ns(~ri, t) µt −(µa+ ∆µa)ns(~ri, t) µt ) = Ii−1·(−)∆µa µt ns(~ri, t)(4.8) So the sensitivity, the Jacobian element, simply reduces to: Ji=−ns(~ri, t) µt (4.9) This, however, does not consider that the sensitivity also depends on the photon’s total path length and how many other voxels a photon passes when going from source to detector as is expressed by the denominator in equation 4.3. We are looking for a way to reduce the sensitivities in voxels that are visited by photons of relatively long path lengths from source to detector. These photons will have small weights when they reach the detector. The sensitivities of those voxels mostly visited by photons of relatively short path length should be larger relative to the rest. There are two ways to account for that. Either we multiply each value dropped by a photon packet in a voxel by the packet’s remaining weight Wwhen it is detected. Or we divide each value by the weight that was dropped on the entire path from source to detector, 1−W. To compare both methods, we plot the ratio of the two different factors for different remaining weights W. The ratio is: Factor Ratio =W 1 1−W =W−W2(4.10) 32 CHAPTER 4. NUMERICAL SOLUTION OF JACOBIAN 0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2 0.25 Remaining Weight W Comparison of Methods Ratio Figure 4.1: Factor of remaining weight divided by factor of division by dropped weight. From figure 4.1 we can see that the choice of method depends on how much weight the photon packets have when they reach the detector. Both methods lead to relatively smaller sensitivity values in voxels far away from source and detector because these voxels are visited by photons that will have relatively low weight at the detector. To what extent the two methods shift the sensitivity ratio between near and far voxels, however, depends on how much weight photon packets still have when detected. If most photons have weights of less than 0.5 at the detector, the method of multiplying by the remaining weight leads to a higher difference in sensitivity values between near and far voxels than dividing by the dropped weight. The effect is reversed if photon packets arrive at the detector predominantly with weights larger than 0.5. The amount of weight photon packets have at the detector depends on the geometry, the optical properties of the materials in the geometry and the source-detector separation distance. In the configurations I worked with, most photons arrive with a weight of less than 0.5 at the detector and the method of multiplying the deposited values by the remaining weight gives better results. The Jacobian element Jiin voxel iis then computed as: Ji=−1 Np Ndp X j Ns,j,i ·1 µt,i ·rWj(4.11) All detected photons Ndp have a different remaining weight rWjand a different 4.2. IMPLEMENTATION 33 number of scattering events Ns,j,i of photon jin voxel i. The total attenuation coefficient µt,i can be different in each voxel i. The division by the total photon number Npserves as a normalization. 4.2 Implementation In the previous section, I explained the theoretical approach to computing the Jacobian within the Monte Carlo method. Here, I want to go into more detail about how to actually implement it in the code. I acquired all the relevant programming knowledge and skills from the excellent book ”A complete guide to programming in C++” by Ulla Kirch-Prinz and Peter Prinz [14]. The central quantity we are interested in is the scattering density ns(~ri, t)of all detected photons. To extract it, we need to record all scattering events. Computationally, the challenge is that we do not know whether a photon packet will eventually be detected. So a photon’s scattering events have to be recorded even though it might not be detected. If the photon packet is terminated without reaching the detector, the information about the scattering events is of no use and hence deleted. For detected photons, the location of all scattering events and time of detection has to be added to the right location of the Jacobian matrix. It is critical but non-trivial to develop a computationally efficient implementation for this process. The code of MBioICFO is written in the object oriented programming language C++, in which each physical object like the source or the detector is defined in a class with its own specific properties. To dynamically allocate memory space to a matrix containing the location of scattering events, I implemented a class to create sparse matrices. The advantage of sparse matrices is that they only require as much memory as is really needed. So when a photon is launched, no memory is yet occupied by the sparse matrix. Only for each scattering event, the memory space is expanded as the coordinates and the total attenuation coefficient in that voxel are saved. On detection or termination of the photon, the sparse matrices occupy differently large memory space according to how many scattering events the photon experienced. Figure 4.2 provides a visualization of this process. At each scattering event, the value 1/µtin that voxel is written to the sparse matrix along with the voxel coor- 34 CHAPTER 4. NUMERICAL SOLUTION OF JACOBIAN Photon launch Information of scattering events is added continuously to sparse matrix x1 y1 z1 1/μt x2 Photon termination without detection Photon detection Memory space of sparse matrix is freed The values in sparse matrix are transferred to Jacobian matrix. The remaining photon weight is considered Sparse matrix Jacobian matrix More photons? Memory space of sparse matrix is freed Yes No All values in Jacobian matrix are divided by total number of photons launched Normalization Output Jacobian x1 y1 Vijk Figure 4.2: Construction of Jacobian in Monte Carlo method. 4.2. IMPLEMENTATION 35 dinates. At some point, the photon is either detected or terminated without detection. In the latter case, the memory space of the sparse matrix is freed and the next photon can be recorded. At detection, the information in the sparse matrix is added to the right place in the Jacobian matrix where the remaining weight of the detected photon is considered. The Jacobian matrix has one element for each voxel and each time bin defined in the input, so its size corresponds to NxNyNzNt, with Nibeing the number of voxels in direction iand Ntthe number of time bins. Note that this size can be specified in the input, so it usually does not contain the full geometry. One can exactly specify where in the geometry and for what time bins the Jacobian is to be constructed. The voxel coordinate is read out from the sparse matrix and the time bin is assigned according to when the photon was detected. The sparse matrix’s memory space is then freed and ready for the next photon. Once the program has simulated all photons, the Jacobian is normalized by dividing all its elements by the total number of photons launched. The user can then use Matlab to read out the separate file reserved for the Jacobian. 42 CHAPTER 5. EVALUATION z [mm] y [mm] Analytical Solution 2.5 5 7.5 10 12.5 5 10 15 20 25 30 −30 −29 −28 −27 −26 −25 −24 −23 −22 −21 −20 −19 Figure 5.7: Analytical solution plotted in logarithmic scale. z [mm] y [mm] Monte Carlo Solution 2.5 5 7.5 10 12.5 5 10 15 20 25 30 −30 −29 −28 −27 −26 −25 −24 −23 −22 −21 −20 −19 Figure 5.8: Monte Carlo solution plotted in logarithmic scale. z [mm] y [mm] Difference MC and Analytical Solution 2.5 5 7.5 10 12.5 5 10 15 20 25 30 −30 −29 −28 −27 −26 −25 −24 −23 −22 −21 −20 −19 Figure 5.9: Absolute difference plotted in logarithmic scale. z [mm] y [mm] Relative Difference 2.5 5 7.5 10 12.5 5 10 15 20 25 30 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 5.10: Relative difference. 5.2. THE JACOBIAN IN THE SEMI-INFINITE MEDIUM 43 The plot of the relative difference between the Monte Carlo and the analytical solution 5.10 suggests good agreement of the two. As for the infinite medium, the relative difference is only large for regions sparsely visited by photons. Note that we are directly comparing Monte Carlo and analytical solution without taking any normalization measures. For a more detailed comparison, we plot the Jacobian along the symmetry axis between source and detector as indicated by the dashed line in figure 5.7. Figure 5.11 shows the unnormalized Jacobian along the symmetry axis averaged over 1 cm in x-direction. Additionally, figure 5.12 shows the natural logarithm of the Jacobian along the same line, also not normalized. So even without normalization measures, absolute values and decay of the Jacobian show reasonable agreement. We can reduce the statistical noise fluctuations in the Monte Carlo solution by simulating an even higher number of photons. For this present study, however, limited computational power did not allow that. 0 2 4 6 8 10 12 14 0 1 2 3 4 5 x 10−8 z [mm] Jacobian Jacobian along symmetry axis between source and detector Theory Monte Carlo Figure 5.11: Comparison between the numerical and the analytic solution. 44 CHAPTER 5. EVALUATION 0 2 4 6 8 10 12 14 −23 −22 −21 −20 −19 −18 −17 −16 z [mm] ln(Jacobian) Natural logarithm of Jacobian along symmetry axis between source and detector Theory Monte Carlo Figure 5.12: Normalized logarithm of numerical and analytic solution. 5.3 The Jacobian at Different Detection Times As mentioned earlier, the implementation developed in the course of this work also allows to find the Jacobian for different time bins of detection. Depending on when a photon is detected, the sensitivities in different regions of the material vary. So looking at the detector signal at different times, we obtain different depth sensitivities. The following figures show the Jacobian in five time bins of the simulation in the semi-infinite medium as presented above. In general, we can freely specify the maximum time and number of time bins in the input. The time bins shown in the figures represent the time when most photons are detected. Clearly, the figures show that for later detection times, the sensitivity is much higher at larger depths and vice versa for early detection times. All five figures are plotted with the same colorscale as depicted by the colorbar in figure 5.13. The plots again show the natural logarithm of the Jacobian averaged over 1 cm in x-direction. 5.3. THE JACOBIAN AT DIFFERENT DETECTION TIMES 45 z [mm] y [mm] 100ps−200ps 2.5 5 7.5 10 12.5 5 10 15 20 25 30 −30 −28 −26 −24 −22 −20 Figure 5.13: First time bin plotted in logarithmic scale. z [mm] y [mm] 200ps−300ps 2.5 5 7.5 10 12.5 5 10 15 20 25 30 Figure 5.14: Second time bin plotted in logarithmic scale. z [mm] y [mm] 300ps−400ps 2.5 5 7.5 10 12.5 5 10 15 20 25 30 Figure 5.15: Third time bin plotted in logarithmic scale. z [mm] y [mm] 400ps−500ps 2.5 5 7.5 10 12.5 5 10 15 20 25 30 Figure 5.16: Fourth time bin plotted in logarithmic scale. 46 CHAPTER 5. EVALUATION z [mm] y [mm] 500ps−600ps 2.5 5 7.5 10 12.5 5 10 15 20 25 30 Figure 5.17: Fifth time bin plotted in logarithmic scale. 5.4 The Jacobian in MRI of Human Head Figure 5.18 shows an anatomical magnetic resonance image (MRI) of a human head 1. The different tissue types are the skull, the cerebrospinal fluid (CSF) and the gray/white matter. To demonstrate the utility of the method developed, we compute the sensitivities at different detection times in this MRI geometry. This shows a potential application of the work presented. In diffuse optical tomography (DOT), several source-detector pairs are placed on the head to measure physiologically relevant variations of the optical properties in the brain [7]. Traditionally, the challenge has been to get enough depth sensitivity to reach the gray and white matter tissue of interest in the brain. For the simulation illustrated in the figures below, we choose a relatively large source-detector separation distance of 3.4 cm to achieve a larger depth sensitivity. The collimated source on the head surface that would be used in the measurement is approximated 1Data courtesy of Yodh lab at University of Pennsylvania 5.4. THE JACOBIAN IN MRI OF HUMAN HEAD 47 Figure 5.18: Image from MRI measurements of a human head. Material 1 is the skull, material 2 the cerebrospinal fluid (CSF) and material 3 and 4 are grey/white matter. by an isotropic source at a depth of approximately 1/µ0 s. The detector is also located at that depth for symmetry reasons. The voxels in the geometry are cubic with an edge length of 1 mm. 109photons were simulated. Figures 5.19 to 5.26 show the Jacobian in the MRI geometry for different detection times. The colorbar in figure 5.19 is valid for the plots of all time bins and gives the natural logarithm of the Jacobian. As pointed out in chapter 3, the Monte Carlo method is important to determine photon transport in media in which the diffusion approximation (see section 2.1) does not hold. It can be seen in figure 5.18 that the diffusion approximation is clearly violated in the CSF, since the absorption coefficient µais almost a fifth of the scattering coefficient µsand the reduced scattering coefficient µ0 sis even a bit smaller than µs. So this is a good example where the diffusion equation 2.2 cannot be used and, accordingly, we have to rely on Monte Carlo solutions. In the figures below, the Jacobian looks very noisy in the CSF. This has mainly two reasons. First of all, as opposed to the solutions presented in the infinite and semi-infinite medium, the Jacobian is shown for only one layer of voxels because the location of tissue types in other layers is obviously different. So we can only take into account photons scattered in this layer. Additionally, the scattering coefficient is very small, two orders of magnitude smaller than in the brain. So there are fewer scattering events contributing to the construction of the Jacobian (see section 4). This leads to more noise. Noise is reduced by a shorter source- 48 CHAPTER 5. EVALUATION detector separation at the cost of smaller depth sensitivity. We can see that there is higher depth sensitivity for larger detection times. While the signal at the detector in the time span from 200 ps to 300 ps is mainly unaffected by the optical properties in the brain, the signal measured for example at 600 ps to 700 ps does contain information of the brain region. We also see how the sensitivities are lower in the 900 ps to 1000 ps time bin, because few photons are detected this late. Figure 5.19: Jacobian in first time bin plotted in logarithmic scale. Sourcedetector separation is 3.4 cm. 109photons were simulated. Figure 5.20: Second time bin plotted in logarithmic scale. Figure 5.21: Third time bin plotted in logarithmic scale. 5.4. THE JACOBIAN IN MRI OF HUMAN HEAD 49 These figures illustrate the value of the code to visualize sensitivities in heterogeneous media at different detection times. The findings from this simulation also agree with a similar study conducted by Boas et al [3]. Figure 5.22: Fourth time bin plotted in logarithmic scale. Figure 5.23: Fifth time bin plotted in logarithmic scale. Figure 5.24: Sixth time bin plotted in logarithmic scale. Figure 5.25: Seventh time bin plotted in logarithmic scale. 50 CHAPTER 5. EVALUATION Figure 5.26: Eighth time bin plotted in logarithmic scale. Chapter 6 Conclusion In this work, I have presented an extension of the Monte Carlo method for diffuse optics. The existing Monte Carlo method and its importance was thoroughly introduced. I explained the significance of finding the sensitivity matrix, also called Jacobian, for any arbitrary geometry to solve the inverse problem and determine a structural image of the tissue. An approach to numerically determine the Jacobian for time-resolved spectroscopy (TRS) in one simulation run was laid out. The method was verified in the infinite as well as in the semi-infinite medium. Besides, I have shown the capability of the program to generate the sensitivities corresponding to different detection times. The implementation has proven to be computationally efficient and the results show good agreement with the analytical solutions for the number of photon packets simulated. I have further demonstrated the power of using the new implementation of MBioICFO developed in this work to determine the sensitivity matrix in any heterogeneous medium at different detection times. Data from MRI measurements can be read in and used as the geometry. The advancement of the Monte Carlo method can stimulate and enhance the power of diffuse optics technology in general. The numerical simulations provide a means to learn more about photon transport in tissue and as the implementations become more computationally efficient, the method becomes even more useful. Hopefully, the Monte Carlo method contributes to making diffuse optics technology a valuable tool in the hospital. 51