Fluorescence lifetime measurements in heterogeneous scattering medium

Abstract. Fluorescence lifetime in heterogeneous multiple light scattering systems is analyzed by an algorithm without solving the diffusion or radiative transfer equations. The algorithm assumes that the optical properties of medium are constant in the excitation and emission wavelength regions. If the assumption is correct and the fluorophore is a single species, the fluorescence lifetime can be determined by a set of measurements of temporal point-spread function of the excitation light and fluorescence at two different concentrations of the fluorophore. This method is not dependent on the heterogeneity of the optical properties of the medium as well as the geometry of the excitation–detection on an arbitrary shape of the sample. The algorithm was validated by an indocyanine green fluorescence in phantom measurements and demonstrated by an in vivo measurement.


Introduction
The fluorescence of the organic molecules, proteins, and nanoparticles is usually sensitive to environment through nonradiative decay process. 1,2 Many applications in biology are using fluorescence probes, which are specifically designed or have intrinsic characters to change their fluorescence properties, in most case through the nonradiative decay pathways, by the environmental changes. [3][4][5] In some cases, the autofluorescence of biological sample reports the environmental change through such mechanism and is very useful in the actual clinical purpose. [6][7][8] The fluorescence intensity measurements have been widely used in bio-and medical sciences. However, the intensity-based fluorescence techniques do not yield absolute quantities and always require a reference to quantify the fluorescence. Furthermore, it is not easy to apply for the scattering system because of the nonlinearity of intensity. On the other hand, lifetime measurements either by the time-domain or the frequency-domain technique is an absolute method to measure the fluorescence decay function and can sensitively monitor the fluorophore environment change. This technique is more robust because the artifact due to the contamination, which is generally difficult to eliminate, is more easily distinguished on the time axis or the phase/amplitude space and many other problems in the intensity measurements do not directly contaminate the data. The actual biological system is highly scattered and significantly heterogeneous. Therefore, in principle, the quantification of the fluorescence lifetime is not simple because the fluorescence temporal response also depends on the time-of-flight in scattering system. 9 The quantification method of the fluorescence lifetime in a homogeneous scattering system was first demonstrated by Patterson and Pogue 10 and widely used in the lifetime estimation for biological applications. However, the actual problem is an inverse problem of the fluorescence diffusion or transport equation for heterogeneous system and thus not just a problem to determine the lifetime; the problem of fluorescence diffuse optical tomography (FDOT) should be solved in principle. 11 This task is very hard and cumbersome. To overcome the mentioned challenges, many simplified models, e.g., using a small distance from the excitation position for the fluorescence measurements to neglect the broadening by the scattering, 12 analytical solution of simple geometry, 13 Green functions of the excitation light and emission light with a simple geometry under Born approximation, 14 or the asymptotic behavior of fluorescence temporal profile, 15,16 have been applied for the extraction of true fluorescence lifetimes. However, the measurements of fluorescence lifetime in heterogeneous scattering medium are a subject of debate and research.
In this contribution, we have proposed another simplified method to extract the fluorescence lifetime function without the actual modeling of the system. The basic concept of the proposed method is that, when the optical properties do not change in the excitation and emission wavelength regions, a very simple equality can be derived between the temporal point-spread function (TPSF) of the excitation light with and without fluorophore interested and the TPSF of fluorescence light. 17 The FDOT using this concept has been known as "total-light algorithm." In this FDOT, the excitation TPSF without the fluorophore is constructed using this equality with the excitation and fluorescence TPSFs and preknowledge on the fluorescence lifetime function. 17 Conversely, the equality means that the fluorescence lifetime function can be determined using the TPSFs of excitation and fluorescence with and without the fluorophore. These concepts and Monte Carlo validation were partly reported in the previous conference. 18 In this paper, we aim to propose a method to determine the fluorescence lifetime using this equality. First, we generalize the theory using the radiative transport equation (RTE), and then validate the theory with phantom experiments. Finally, we demonstrate this method with an in vivo measurement.

Theory
Let us consider a fluorophore object in a scattering medium.
The main concept of the analysis is based on the energy transaction from the excitation light to fluorescence light. This concept was originally described in a previous paper for homogeneous systems 10 and here we will follow this but extend to heterogeneous systems. The fluorescence is generated by the absorption, with the energy conversion ratio γ, which is so-called the quantum efficiency. Once the excited state is generated, the fluorescence photon is emitted with a delay due to the statistical nature of the excited-state kinetics. The excitation light simultaneously lost its energy at the excitation process by the fluorophore absorption. If there is an object that contains the fluorophore, the only arriving photon to path though the region is contributing for the absorption-emission process. The excitation-fluorescence conversion causes the wavelength and the directional distribution changes; the absorption process has a memory of the direction, i.e., maintaining the original direction, but the fluorescence process works an isotropic directional change when the rotational relaxation of the fluorescence molecule is sufficiently fast; more accurately, many fluorescence process are involved during the measurement and statistically homogenized the rotational relaxation of the fluorescence molecules.
If the optical properties of the medium do not depend on the wavelengths at the excitation and emission, the trajectories of the photons through the object region will not change under the validity of the diffusion approximation. Consequently, if the fluorescence decay is infinitely fast, the fluorescence TPSF of this system, the instantaneous fluorescence TPSF (IF-TPSF), will be same to the TPSF of the excitation light, which has been absorbed by the fluorophores. Therefore, the change of the TPSF of the excitation light with and without the fluorophores is the exact same as the TPSF of absorbed excitation light, which can be estimated from IF-TPSF with knowledge of quantum efficiency. Finally, an equality between the change of TPSFs of the excitation light and IF-TPSF is derived. This idea can be proven from RTEs.

Fluorescence Radiative Transfer Equation and
the Concept of the Analysis In this section, we describe the time-dependent fluorescence photon transport in a multiple scattering system. The fluorescence transport equation is a coupled equation depending on the excitation photon transport. It can be decoupled, and the problem becomes two individual transport problems under some assumptions. This statement was proven in case of diffusion equation (DE) system. 17 This fact can be derived from RTEs with the homogeneity of the fluorescence decay profile F: fluorophores can distribute inhomogeneously in space but in a homogeneous local environment around the fluorophores, showing an identical fluorescence decay kinetics. First, we consider that the sample is illuminated by the excitation light Sðr; s; tÞ, where r; s, and t are the position, direction, and the time, respectively. When the most general expression of the photon transport in a strong turbid medium is RTE and the radiance Lðr; s; tÞ can be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 7 5 2 1 c ∂ ∂t þ s · ∇ þ μ a ðrÞ þ μ s ðrÞ Lðr; s; tÞ Pðs; s 0 ÞLðr; s 0 ; tÞdΩ 0 þ Sðr; s; tÞ; (1) where μ s ; μ a , and c are the scattering and absorption coefficients, and speed of light, respectively. Pðs; s 0 Þ is the scattering phase function, which is the probability of the scattering with an incident direction s 0 and the scattering direction s. The right-side integral with respect to the solid angle dΩ 0 sums up all incoming light at this point. In case of the fluorescence, the source term S becomes E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 6 1 2 where FðtÞ, γ; ε, and NðrÞ are the fluorescence decay function, quantum efficiency, molar absorption coefficient, and concentration of the fluorophore, respectively. The suffix "ex" means the excitation light while the "fl" means fluorescence. This source term is neglecting a reabsorption-emission process, where fluorophores are excited by the fluorescence. The convolution integral can be eliminated by introducing L Ã fl , which sat- where μ a;fl ¼ eNðrÞ, meaning the absorption coefficient of the fluorophore. The absorption coefficient μ ex a ðrÞ is separated by two parts, those of the fluorophore μ ex a;fl ðrÞ and medium μ ex a;med ðrÞ. Then, a quantity L v ¼ L ex þ 1∕γL Ã fl at r, s, and t, is introduced. L v satisfies an RTE E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 2 9 4 1 c when μ fl a ¼ μ ex a;med , μ fl s ¼ μ ex s , and P fl ðs; s 0 Þ ¼ P ex ðs; s 0 Þ. These conditions mean the equality in the optical properties at the excitation and fluorescence wavelength. The last term of Eq. (4) comes from the memory of the direction of the excitation light; the direction of the excitation light does not change and only the energy is decreased when the light is absorbed by the fluorophores. The excitation light through the object region works as the source of the light at the point with no fluorophore and maintains the direction. On the other hand, the fluorescence light, generated here, is emitted in isotropic directions when the rotational relaxation of the fluorescence molecules is sufficiently fast to randomize the direction of the fluorophore dipole. Eventually, this difference between the excitation light and fluorescence causes the last term. If the fluorescence target is sufficiently far from the detection point, the difference due to the directional memory is lost and thus this term will be dropped. Assuming can be yielded. The detail derivation can be found in Appendix. Under the DE, the directional memory of the light is completely ignored and this relationship is correct with the independency of the optical properties on the wavelength. The derivation here is a more generalized version of the previous paper. 17 It is worth noting that this equality is always satisfied in any case under these assumptions discussed above. The assumptions do not include the inhomogeneity of the medium, geometry, and boundary conditions. In addition, this result is also valid in case of the frequency domain. The FDOT application with this equality has been demonstrated. 19 This method is actually DOT to visualize "absorption" of the fluorophores and, thus, the method can be named as "fluorescence-assisted DOT" (FA-DOT).

Fluorescence Decay Analysis
When the excitation TPSF with and without fluorophores, L ex ðtÞ and L 0 ex ðtÞ can be measured, L Ã fl ðtÞ is yielded by L 0 ex ðtÞ − L ex ðtÞ from Eq. (6). Therefore, referring to the previous definition of IF-TPSF L Ã fl , a simple relationship can be obtained as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 3 6 0 L fl ðtÞ ¼ γ In time-domain measurements, all measured temporal responses are broadened by the instrumental function (IRF) H, which is determined by the system. The actual measurement temporal profiles I fl ðtÞ, I 0 ex ðtÞ, and I ex ðtÞ can be expressed by convolutions L fl ðtÞ Ã H, L 0 ex ðtÞ Ã H, and L ex ðtÞ Ã H, respectively, when IRFs are the same in these measurements, where "*" means the convolution operator. Commuting the convolution operators, one can finally yield a simple relationship E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 6 3 ; 2 2 9 I fl ðtÞ ¼ α½I 0 ex ðtÞ − I ex ðtÞ Ã FðtÞ; where α is a proportional constant determined by the quantum yield, difference of the sensitivity detector and other differences, such as the difference of losses of optics, between the excitation and fluorescence measurements. This equation is very much similar to the conventional fluorescence temporal profile for the clear solution but ½I 0 ex ðtÞ − I ex ðtÞ works as an instrumental response function in conventional decay analysis. Therefore, the common fluorescence lifetime fitting analysis can be applied for Eq. (8); just replace the IRF by ½I 0 ex ðtÞ − I ex ðtÞ in the fitting program. It is worth noting that the equation does not include IRF explicitly and, therefore, one can use the expression with the measurement profiles without the knowledge of IRF.
The relationship can be extended to the measurements for changing the fluorophore concentration. Two measurements with different concentration of fluorophore I a and I b , one can modify Eq. (8) as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 3 2 6 ; 7 0 8 I a fl ðtÞ − I b fl ðtÞ ¼ α½I b ex ðtÞ − I a ex ðtÞ Ã FðtÞ: (9) This expression may particularly be useful in measurements where the zero concentration of the fluorophore cannot be measured and there is a background.

Monte Carlo Simulation
The directional memory term, Eq. (5), is omitted to derive L 0 ex ¼ L ex þ 1∕γL Ã fl from RTE, in contrast to that from DE, which does not need to consider the direction of light. Once this equality is attained, Eqs. (8) and (9) are automatically satisfied with the assumption of homogeneity of the fluorophores, i.e., homogeneous lifetime and quantum efficiency, and the invariance of IRFs. Therefore, the validity of this equality is the key in our fluorescence lifetime analysis. Here, we conducted Monte Carlo (MC) simulations to confirm the equality. In particular, we focused on a geometry with the smallest source-detector distance (null distance), where DE is not applicable.
The simulation was using a general-purpose graphics processing units (GPGPU) with a modified program based on CUDAMCML coded by Lund group. 20 The main differences of the program are following: (1) TPSF is calculated by histogram of the photons according to the path length traveled in the medium; (2) the photon energy only takes a binary values 0 or 1; (3) a sphere fluorescence object is set in an infinitely large semiinfinite medium. The actual simulation geometry is shown in the inset of Fig. 1(a). In the simulation, the unity quantum efficiency and the infinitely short fluorescence lifetime were assumed. This simplification does not affect the validation of the equality. The quantum efficiency γ linearly affects the scale of fluorescence intensity and this intensity change is compensated by the prefactor 1∕γ in the equality. The statistics of MC results becomes worse with a smaller value of γ but does not change the shape of the TPSF. The fluorescence TPSF is given by the convolution with fluorescence lifetime function in case of the finite lifetime. For the validation of the equality, the deconvolution procedure of the fluorescence TPSF with the fluorescence lifetime function is needed, resulting in just an increase of the numerical error. Therefore, the simplification does not affect the validation of the equality.
The optical properties of the medium were similar values of the phantom; μ 0 The optical properties of the object were the same to the medium except the absorption due to the fluorophores 0.039 mm −1 . The scattering anisotropy g was set to 0.0 or 0.9, therefore, μ s ¼ 1 mm −1 or 10 mm −1 to maintain the same value of μ 0 s ¼ 1 mm −1 , respectively. The program was run on a GPGPU (Tesla C1060, C2025 and C2020, Ndivia) and 5 × 10 10 photons were traced.
The TPSFs by the MC simulation with g ¼ 0.9 and with a 6-mm sphere target at (0 mm, 0 mm, 10 mm) are shown in Fig. 1(a). The injection and detection points were both at the origin. The excitation light with the target object decayed more quickly due to the absorption of the target. The fluorescence TPSF was delayed compared with the excitation decay, because the fluorescence occurs after the excitation light reaches the object. The theory predicts the relationship for the intensities In this simulation, the quantum efficiency of the emission was unity γ ¼ 1 and the fluorescence lifetime was ignored. Therefore, L Ã fl is identical to L fl in the figure. The sum L ex þ L fl is shown by L v in the figure, but L 0 ex cannot be visually distinguished from L v . Figure 1(b) shows more clearly the differences between L v and L 0 ex for the simulation with z ¼ 5, z ¼ 10, g ¼ 0, and g ¼ 0.9. No systematic deviation can be found in both results at z ¼ 10. On the other hand, a small deviation can be found in the result at z ¼ 5 with g ¼ 0.9. This can be attributed to the directional memory of the excitation light. This directional memory should be lost for the assumption of the theory since fluorescence does not have such memory. In case of g ¼ 0.9, the directional memory still remains at 5 mm because of insufficient number of scattering events. On the other hand, the memory is lost even with single-scattering event in case of g ¼ 0.0 and, therefore, the difference between excitation and fluorescence paths is quickly reduced with the 5-mm separation. The MC results with different geometries are presented in Appendix and are consistent with the above result. Therefore, the equality is valid when the fluorescence target is sufficiently far from the detection point, i.e., a sufficiently long distance from fluorophores where the diffusion approximation for fluorescence light is valid. It is important to note that RTE is needed to describe excitation light.

Experimental System
The block diagram of the experimental setup, which is similar to a previous one, 19 is shown in Fig. 2. Briefly, a mode-locked laser (Tsunami or Mitai, Spectra Physics) at 780 nm was coupled to a polarization-maintained single-mode fiber (P1-780PM, ThorLabs). The pulse repetition was reduced by an electric optical modulator for Tsunami or an acoustic optical modulator for Mitai. The fiber was bifurcated with a ratio of 1:99 (FC780-99B-FC, ThorLabs) for a power monitor (PM100USB, ThorLabs). The main end of the fiber splitter was connected to a multimode (MM) fiber switch (FOSW-1-8-L-62-L-2, Lightwave Link Inc.) for changing the excitation point. The output of the fiber switch was connected to another MM fiber and the end of the fiber was connected to a collimator (F230FC, ThorLabs) to illuminate a sample. The collimator was not used in phantom measurements and the end of the fiber directly contacted on the surface of plastic phantom. A 3-mm bundled fibers with or without a collection lens was located at appropriate detection position. Here also, the lens was eliminated and the fiber contacted to the phantom in the phantom measurements. One of the bundled fibers was selected by a custom-designed fiber switch and the output was divided into excitation and fluorescence beams by a dichroic mirror and then each beam was filtered to select excitation (FF01-780/12, Semrock) or emission (FF01-834/LP, Semrock) wavelength. Finally, each beam was detected by a multichannel plate photomultiplier tube (R3809U-61, Hamamatsu) connected to a timecorrelated single-photon counting (TCSPC) system (TRSF20F, Hamamatsu). The count rate was measured at the output of a constant fraction discriminator by a fast counter board with a 1-s sampling time. Since the dead time of the TCSPC system was not negligible, the count rate measured by the fast counter board was always higher than the event rate of the TCSPC system. Tsunami and the Mitai laser systems were used for the measurements of phantoms and in vivo, respectively. Approximate injection and detection points on the abdominal region in the in vivo measurements are indicated by red and green circles as shown in the inset superimposed picture of Fig. 2. In the phantom measurement, the excitation and fluorescence were sequentially measured with a single-channel system. The excitation and fluorescence were selected only by narrow-band filters. On the other hand, the excitation and fluorescence were, simultaneously, measured in the animal measurement.
The fluorescence lifetime was analyzed with a multiexponential function by Levenberg-Marquardt method to minimize χ 2 .  diluted from 10% Intralipid solution (Fresenius Kabi AG) by deionized water. The blank solution was an 1% Intralipid solution. An ICG aqueous solution with millimolar concentration was prepared and, after the preparation, quickly diluted about 100 times in two different solutions with 1% Intralipid and water. The water solution was quickly measured after the preparation by an absorption spectrometer to estimate the concentration of the initial ICG solution. We used the molar absorption coefficient 0.16 cm −1 μM −1 of ICG aqueous solution at the peak of the absorption spectrum. Finally, the appropriate dilution by 1% Intralipid solution was made to prepare ICG-Intralipid solutions. Once ICG mixed in Intralipid, the solution became very stable. The measurements were conducted at about the middle plane of the cylinder. The cylinder could be rotated and the origin of the angle was defined at the closest point from the target.

Resin Phantom
The fluorescence lifetime of the ICG-Intralipid complex solution was determined by separate measurements for validation of the proposed method. The ICG-Intralipid solution was measured in a thin glass cell with about 1-mm thickness to eliminate the effect of the multiple scattering and measured by the same TCSPC system. The decay was expressed by a single-exponential function with the lifetime 0.63 AE 0.04 ns in a range from 1 to 5 μM ICG in 1% Intralipid solution. The value was also confirmed at a very low concentration of the mixture solution by a confocal microscope system with a TCSPC system. The lifetime value is in good agreement with the published value. 12

Animal Measurement
Animal measurement was conducted under the regulation of animal ethic rule by Hokkaido University. A male hairless rat (8 weeks of age, body weight about 250 g, HWY/Slc, Japan SLC) was anesthetized by urethane (2.0 g∕kg, i.p.). The animal was held upside down and measured from the abdominal surface as shown by the picture in Fig. 2. The excitation and detection optics were fixed to a fiber holder, touching the skin surface, with a black sponge to reduce contamination from the surface scattering. The excitation and detection points and the middle point of them are indicated by the circle, box, and cross symbols in the picture, respectively. The separation between the excitation and detection points was 12 mm. The animal received a bolus injection of a 0.1-ml of 200 μM ICG solution from the tail vein. ICG flew in blood vessels and accumulated in tissue. After the measurements, the rat was euthanized with an overdose of Nembutal.  Table 1. The TPSFs and the fitting (data not shown) were very similar to the results in Fig. 3 and the lifetime was almost constant except at the lowest concentration (0.2 μM).

Phantom Measurements
The difference of the excitation TPSF with and without the target was too small, and the variation of data points at the lowest concentration might cause the large deviation. These fluorescence lifetimes were close to the value obtained by a separate standard measurement and published values 12 , but both values are smaller than those values. Further, the difference of the value with the geometry Fig. 3(b), where the detection position from the target was longer, is more significant. Actually, the estimated lifetime became shorter with an increase in the distance from the target position (data not shown). This can be explained by a violation of the assumption in the theory. The wavelength dependence on TPSF was measured with a larger POM resin (data not shown) and TPSF was analyzed by fitting to the semi-infinite solution of DE with a hybrid approach of extrapolated boundary condition and partial current boundary condition to estimate the scattering and absorption coefficients. 21 In the analysis, the refractive index was assumed to be 1.5. The difference of the scattering coefficient in the range of the excitation and emission was very small (1.0 mm −1 at 760 nm and 0.98 mm −1 at 820 nm). On the other hand, the absorption in the emission wavelength range was significantly higher than that in the excitation (2.1 × 10 −4 mm −1 at 780 nm and 7.5 × 10 −4 mm −1 at 820 nm). This suggests that the fluorescence TPSF was narrower than expected because the higher absorption exponentially eliminates the longer paths resulting in smaller dispersion of the path length. Moreover, the longer distance from the fluorescence target makes more significant narrowing of the path length distribution. Therefore, the result was qualitatively explained by the difference of the absorption coefficient.
A simple moment analysis can be conducted for the quantitative discussion of this systematic deviation. The first moment of the fluorescence TPSF is given by a sum of the fluorescence lifetime, the mean transit times to/from the target and the mean IRF as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 3 8 9 where the superscripts distinguish between the excitation and fluorescence light and IRF. The suffixes "SO" and "OD" indicate the path from source to object and object to detector, respectively. The fluorescence is assumed to decay with a singleexponential function. If the optical properties are constant in the excitation and fluorescence ranges, ht fl OD i ¼ ht ex OD i is valid. The term on the right-hand side can be calculated from the difference of TPSFs with and without the object. The violation of the constant optical properties in the range of excitation and fluorescence breaks this equality. The difference Δτ fl of the moments E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 3 2 6 ; 2 3 6 can estimate the difference of ht fl i from that expected by the constant optical properties assumption, and, eventually, this value gives the deviation of the estimated fluorescence lifetime from the true value. To calculate the difference, we need to solve explicitly the RTE or DE with the knowledge of the optical properties. In a simple case, like a point fluorescence object and a homogeneous medium with a simple geometry, one can use the analytical solution. Here, an analytical solution of cylindrical geometry under DE 21 was employed and the above optical constants of POM were used in the calculation of TPSFs. The source of the light was set at the center of the object and the first moment of TPSF at the detection point was numerically  calculated from the TPSF. The moments with the geometry in Fig. 3(b) were estimated as ht fl OD i ¼ 1.181 ns with the optical constants at 820 nm and ht ex OD i ¼ 1.264 ns with those at 780 nm. The moment of the fluorescence TPSF is 0.083 ns shorter than the value expected under the assumption of the proposed method. In the same way, the difference is 0.027 ns with the geometry in Fig. 3(a). The differences in experimental results between the standard method 0.63 ns and the average of our proposed method 0.56 ns with Fig. 3(b) and 0.61 ns with Fig. 3(a) (exclude the data at 0.2 μM) are 0.07 and 0.02 ns, respectively, and these values are in very good agreement with the estimation.
In general, when the assumption on the optical properties in the range from the excitation to the emission is not valid, the analysis of the trajectory from the object for both of the excitation and emission is required. Eventually, the problem returns to the FDOT problem. This systematic deviation from the true value may be corrected in particular case by some assumptions.
The over-or underestimation of the fluorescence lifetime is caused by the difference of the absorption and scattering coefficients. This can be minimized if the fluorescence target is closed to the detection point, to shorten the path length of fluorescence.

In Vivo Fluorescence Lifetime Measurements
The relative excitation and fluorescence intensities are shown in Fig. 4. The ICG solution was injected at time 0. The relative intensities were calculated by the ratio to average intensities of about the first 10-min period shown by the filled symbols at 328 s. Here, the time point was defined by the middle time of the measurement period. The intensities after the injection were higher than that before the injection. This is due to an artifact of the positioning and contact condition of fiber probe caused by the procedure of injection. Therefore, the intensities could not be compared accurately to those before the injection. A weak background signal from tissue was observed before the injection. This signal was originated from the autofluorescence and the contamination from the excitation light. The signal was not negligible at the later time points. The fluorescence intensity was increasing after the injection, indicating the accumulation of ICG in the tissue. In contrast, the intensity of the excitation light did not change significantly, but was slightly decreasing after the injection due to the increase of absorption by ICG. The time scale of the increase of ICG fluorescence was slow in this experiment. This is probably because the accumulation of ICG in the organs, liver, kidney, or bladder was measured.
The amplitude of the TPSFs was calibrated before the comparison with the reference profile, which was the first measurement data after the injection. The calibration was conducted as follows. First, the total event number was normalized by the count rate measured by the fast counter to calibrate the dead time loss of the TCSPC system and the difference of the actual accumulation time. Then, the amplitude of the profile was corrected by the monitored power to compensate the excitation power fluctuations. Finally, the area of each profile was scaled by a scale factor of the reference to make the same amplitude axis to the reference one. This procedure is very important because the excitation profile change is usually very small and a small error in the calibration may cause a large error in the comparison.
The temporal profiles are shown in Fig. 5. The changes of the excitation TPSFs cannot be visualized in the semilog plot. However, the relative change from an initial TPSF at 328 s. after the injection can actually show a gradual decrease of the profile with time as shown in the inset of Fig. 5(a). In principle, the difference TPSF from the profile before the injection can be used, but the first TPSF data after the injection was employed because of the artifact above explained. In contrast, the fluorescence TPSFs simply show a gradual increase with time and become a more monotonic decay curve as shown in Fig. 5(b). The initial profile before the injection indicates that a background existed and all profiles have the contamination from this background signal. On the other hand, the difference profiles of the fluorescence TPSFs show very similar slopes and gradually increase with time as shown in the inset of Fig. 5(b). This suggests that the background signal was constant during the measurements and the difference of the profile indicated the accumulation of ICG. Figure 6 shows the data at 2616 s and the best fit to Eq. (9) with a single-exponential decay function. From this, the fluorescence lifetime, 0.63 ns, could be estimated. The fitting results are depending on the statistics of data. Therefore, the intensity scale was further calibrated by the accumulation time of measurement to reflect the actual statistics of data. The fitting curve seems to have a bit faster slope than the experimental one. This is most probably due to the problem of the fitting weight. Usually, TCSPC data are analyzed by the least-squares fit, assuming a Poissonian variance for the weight of the evaluation function. This usually causes a systematic deviation at a lowcount region. 22 Further, the subtraction data of two Poissonian TCSPC data for the fitting no longer follow the Poissonian statistics. Therefore, the fitting was biased by the incorrect weight of data, particularly at a low-count region. Therefore, more accurate procedure to manage the statistics is needed to improve the fitting. The average fluorescence lifetimes estimated at four different times after the injection was 0.60 AE 0.05 ns. This fluorescence lifetime was very close to a short lifetime component of 0.56 ns of ICG fluorescence decay in plasma.

Discussion
The fluorescence TPSF is dependent on both the fluorescence decay function and the path length distribution of the excitation and fluorescence light. The path length is depending on the system, such as the shape, heterogeneity inside, and measurement geometry, resulting in the estimation of the fluorescence lifetime difficult in general. In this paper, we have proposed a simplified method using changes of TPSF of the excitation and fluorescence light by a concentration change of fluorescence chromophores. This idea has come up from FA-DOT. We have started an RTE formalism rather than DE to derive the simple equality as Eq. (6) with the additional assumption Eq. (5). This comes from the difference of the directional memory of the excitation light and fluorescence. The directional vector of the excitation light from the object should be randomized to satisfy the assumption. This is one of the requirement conditions in the derivation of DE from RTE. However, this does not expect that the excitation light from the injection point need to be satisfy the condition of DE. The MC simulation shows the result with null separation between the source and detector and the results indicate that the equality is maintained except at the shallow object depth. Therefore, the scheme presented here is not completely identical to DE. The proposed method is basically not depending on the heterogeneity nor the geometry of the sample with the assumption of the same optical properties in between the excitation and emission wavelength regions. Under this assumption, the trajectory of the excitation light eliminated by the fluorophore absorption is same as that of the fluorescence. Therefore, the method presented here has an advantage in the application for biological systems, which are highly heterogeneous and has an arbitrary shape. This method does not require any reconstruction nor modeling. Of course, the assumption is presumably not good in general in biological systems, and the failure of the assumption causes a systematic bias in the results as shown in POM experiments. This is not only the problem in time-domain measurements but also in other measurements, e.g., the frequency-domain measurements. This effect can be reduced using a detection close to the fluorophores. It is worth noting that many works on the fluorescence do not consider the differences of the optical properties in the excitation-fluorescence range, meaning that they have same problems here.
The similar advantage is the measurements with the background signal, which are from autofluorescence and contamination from the excitation. The background without the target was negligible in the phantom measurement. However there are some cases in which the background signal is not negligible as shown in the in vivo measurement. The contamination from the excitation can be reduced by the better rejection filters at the detector. The autofluorescence is actually reduced in a longer wavelength range of near-infrared region and might be reduced by manipulation of ingredients in foods for in vivo experiments. However, the deep tissue fluorescence measurements require a very sensitive measurement technique, such as the photon counting method, and it is really difficult to completely eliminate the background from the autofluorescence from massive tissue surrounding the fluorescence target and the contamination from the excitation. The autofluorescence is usually attributed to some chromophores in tissue and will be constant in a time scale of the measurements. The contamination from the excitation will also be constant. Thus, the method here can automatically remove the overlapped background signal because the change of the TRSF is used. One may analyze the fluorescence lifetime by using change of the fluorescence with a solution of DE or RTE, but an explicit modeling is required.
The disadvantage of this method is the measurement on the small change of the excitation light. The fluorescence method is usually very sensitive but this method uses the change of the excitation through an absorption change by the fluorophores. Eventually, this method requires a relatively higher amount of change of the fluorophores to make the change of excitation light visible. The results will be extremely sensitive to the fluctuations of measurements because of the small change of the excitation. The fluctuation of the background signal is also an important factor. In the in vivo measurements, the TPSF change was about 10% at maximum. The blood volume of rat is about 16 ml 23 and the ICG concentration was about 1.25 μM in blood vessels if the ICG solution was homogeneously diluted. But, the effective concentration of ICG was in the order of 0.1 μM estimated from the body weight. The gradual fluorescence increase suggests the accumulation of ICG in tissue rather than the initial circulation of ICG after injection. This suggests that the 0.1-μM order of ICG was below the sensitivity of detection and probably the lowest detectable concentration was an order of μM.
The analysis with Eq. (9) was directly conducted with a single-exponential fluorescence decay function by Levenberg-Marquardt algorithm to minimized χ 2 , which assumes a Gaussian error statistics of data. However, the TCSPC data itself is considered to be Poissonian and further the Poissonian statistics is not valid for the difference of two TPSFs. This causes the systematic error of fitting. More theoretical work for the improved algorithm is required in future works. The other method, Phasor analysis, may also be useful because Eq. (9) is identical to common fluorescence decay fitting functions. 24 Phasor analysis is just the Fourier-domain representation of the temporal response function but very intuitive.

Conclusion
We have demonstrated a fluorescence lifetime measurement method in multiple scattering system. Our primary target application is a medical diagnosis through the fluorescence lifetime change in tissue for cancer detection. However, this method is not limited to such biological system but also applicable to the monitoring of the fluorescence in turbid systems, such as colloid solutions. The fluorescence lifetime measurements are usually very difficult in strong scattering media and this method shall be a simple option of the measurement.
Appendix A: The Derivation of the Equality Eq. (6) The excitation light L ex is satisfied by Eq. (2) and more explicitly expressed as The absorption term μ ex a can be divided by two parts μ ex a;med and μ ex a;fl , where the absorption coefficient of medium and that of the object (fluorophore) at the excitation wavelength. Then, let us introduce a quantity L v ¼ L ex þ 1∕γL Ã fl . Replacing L ex by L v − 1∕γL Ã fl in Eq. (12) except the term related to μ ex a;fl yields an equation as is obtained and then using an integral form of μ ex a;fl L ex with respect to the direction, Eq. (4) can be derived. Further, using the assumption of Eq. (5), one can obtain This equation is the exact same to the RTE for the excitation light in the medium without the object. Therefore, fl . The first term of the assumption of Eq. (5) means that the generation of the fluorescence is isotropic. On the other hand, the second term of the equation means that the excitation light does not change the direction.

Appendix B: MC Simulation with Different Geometry
The MC result with a different geometry is presented in Fig. 7. The photon was injected at the position (−1 cm, 0 cm, 0 cm). Other simulation parameters were same to the previous simulation. The TPSFs of the excitation with and without the target and that of the fluorescence at y = 0 are shown by pseudocolor images. The ratio of L v ¼ L ex þ L fl to L 0 ex is also shown. The white region of the image is the region where the photons do not come due to the finite velocity. A small number of photons came to the region close to the white region as shown by the dark blue color. The red-color region of the ratio coincides this region, but this is caused by a technical problem of the color image. The actual distribution of the ratio distributed almost equally around 1 but the fluctuation is high due to the low statistics. Therefore, the ratio is distributing around 1 in time and space, indicating the equality as Eq. (6) is maintained. The target distance from the source is longer than the previous simulation. This suggests that the directional memory of the excitation light loses before reaching the target and the directional memory term in Eq. (4) can be omitted.