Abstract
We describe a new light transport model, which was applied to three-dimensional lifetime imaging of Förster resonance energy transfer in mice in vivo. The model is an approximation to the radiative transfer equation and combines light diffusion and ray optics. This approximation is well adopted to wide-field time-gated intensity-based data acquisition. Reconstructed image data are presented and compared with results obtained by using the telegraph equation approximation. The new approach provides improved recovery of absorption and scattering parameters while returning similar values for the fluorescence parameters.
© 2011 Optical Society of America
1. Introduction
We recently reported [1] the demonstration of tomographic fluorescence lifetime imaging (FLIM) in vivo to read out of Förster resonance energy transfer (FRET) [2, 3] in live mice, combining diffuse optical tomography [4, 5] with FLIM [6, 7, 8, 9, 10, 11]. By implementing wide-field time-gated image data acquisition of excitation radiation and fluorescence excited by ultrashort laser pulses in a tomographic imaging system [12, 13, 14, 15, 16] and applying an image reconstruction algorithm based on the telegraph equation (TE), we could recover the three-dimensional (3D) spatial distribution of fluorescence quantum yield and lifetime throughout a specimen. In this work, we continue our development of 3D FLIM tomography for reading out FRET of biomolecular interactions in live mice by introducing a new technique for tomographic image reconstruction of weakly scattering samples that utilizes a novel, computationally inexpensive, approximation to the radiative transfer equation (RTE) that replaces the TE. The performance of these two light transport models is compared for experimental data acquired from fluorescence proteins expressed in mouse legs [1]. Briefly, the hind leg anterior tibialis muscles of live mice were transfected by electroporation, and in vivo time-resolved images of the donor [enhanced green fluorescence protein (eGFP)] fluorescence were acquired using a wide-field time-gated imaging system incorporating a gated image intensifier. The mouse legs were transfected with either free donor expressed in the presence of free acceptor (mCherry) or with a tandem FRET construct (eGFP–mCherry) in which the donor was directly linked to acceptor. While the use of the TE successfully distinguished the FRET and non-FRET samples, the reconstruction of the optical properties was less successful.
The TE is a variant of the diffusion approximation (DA), which is a statistical description of photon transport that does not consider the propagation of light according to geometric optics since it assumes that all photons are scattered from their original ballistic trajectories. Unfortunately, the DA is not a good approximation for small highly scattering objects whose size is comparable to the photon’s mean free path length and, therefore, from which photons have a finite probability to leave without any scattering or absorption events. We believe that realistic light transport models for such small scattering objects should take account of both photon diffusion and rectilinear light propagation. These properties are described by the RTE [17], but the numerical solution of the RTE is highly challenging, with prohibitive computational overheads for most applications. In this paper, we suggest a new, computationally inexpensive approximation to the RTE that combines elements of light diffusion and ray optics.
Including light propagation along rays into the light transport model is important for interpretation of images taken by a CCD camera. Images are formed from the energy of absorbed photons by pixels of the CCD array within the camera’s exposure time. Photons reach the CCD array by propagating along rays leaving the surface of a scattering object. Thus, the intensity of light rays per pixel area and per exposure time gives the absorbed energy by a pixel of the CCD array. Absorbed energy is further transformed into an image. Therefore, the time-gating technique relies on intensity-based measurements. The DA, for instance, describes the energy density of radiation, which is proportional to the first moment of the intensity [18]. Previously, we utilized the method of Schwarzschild–Schuster [17] for an image’s transformation from intensity measurements to energy density measurements [19]. However, this method has its limitations. A model describing the intensity of light rather than its moments has a much broader range of applicability.
The image reconstruction algorithm is derived by employing a variational framework in the Fourier domain. It combines the backprojection algorithm, which is a well-adopted reconstruction technique in computerized tomography, and diffusive inverse imaging. It is well suited for wide-field time-gated intensity-based detection. Cartesian meshes are constructed from shadowgrams taken of mouse legs with a voxel edge matching the CCD camera’s pixel size.
This paper is organized as follows. Mathematical and computational details are described in Section 2, where Subsection 2A is devoted to the direct problem of light transport, while Subsection 2B describes the image reconstruction approach. Experimental details and the mesh construction procedure are outlined in Section 3. In Section 4, we compare reconstruction results based on the suggested model with results obtained by using the TE.
2. Methodology
2A. Direct Problem
Light transport in scattering media is modeled by the RTE [17]. Although the RTE provides an almost complete and very general description of the light scattering phenomenon, from a practical point of view, it is too complex and numerically expensive. Computational overheads of numerical solution of the RTE are undesirable in the context of the inverse problem. In this subsection, we consider an inexpensive approximation to the RTE, which is applicable to small highly scattered objects. Our method is based on the TE approximation and directly takes into account the probability of the quantum exit of a photon from the medium in a given direction [20, 21].
To avoid introduction of additional unknown parameters, we assume an isotropic scattering in the medium. This can be justified by introducing an effective scattering coefficient by making the scattering photon’s mean path longer. Then, the RTE in the Fourier domain (ω) can be written in the form of the first-order partial differential equation:
In Eq. (1), I is the intensity of light, is the attenuation coefficient, which is given by , where μ is the transport coefficient, and c is the speed of light. The transport coefficient is reciprocal to the mean path length of a photon in the medium and represents the sum of the scattering coefficient, , and the absorption coefficient, . The vector s is the unit vector in a given direction. In the source term, Eq. (1), u denotes the average intensity:In order to treat Eq. (1) as a first-order partial differential equation, the average integrated intensity, u, must be supplied self-consistently with the RTE. In highly scattering media, this quantity is computed by solving the TE [18], which results from the first two moments of the RTE. That is, multiplication of Eq. (1) by 1 and by s and consequent integration over the whole solid angle leads to the Helmholtz equation, which is the Fourier image of the TE:
where the Helmholtz differential operator is given by where the diffusion coefficient κ is expressed in terms of the attenuation coefficient as In the source term, represents the average intensity of the direct excitation radiation , which satisfies the equation where is the amplitude and is the location of a laser source, which injects excitation photons into the medium in the direction . Thus, the direct excitation radiation is given by which represents an attenuated laser’s ray entering the scattering domain.Because the variation of the refractive index in the medium is neglected, the source location, , and the observation point, r, are connected by a line . Computationally, it is beneficial to compute the intensity of radiation leaving the scattering object by integrating from the observation point. Therefore, the general solution of the RTE [Eq. (1)] at the observation point r is provided in the following form:
where denotes the maximum distance of a ray’s path contributing to the intensity. The exponent in Eq. (7) is proportional to the probability for a photon to exit the medium along s without scattering and absorption events. This probability is integrated with the source distribution along a ray’s path, resulting in the observed intensity. In our experiments, no ballistic photons were observed and, therefore, terms containing are neglected in Eq. (7).Fluorescent light transport is analogous to excitation. Usually, fluorescence is considered as a re- emission of the excitation radiation, which was absorbed by a fluorophore. However, in accordance with our approximation, it can be viewed as the resonantly scattered excitation radiation. The resonant scattering by fluorescent particles is isotropic and accompanied by Stokes shift. The transport of such scattered excitation radiation is governed by the equation
where In Eq. (9), η denotes the quantum yield and τ is the lifetime, which together characterize a fluorophore’s response to the excitation radiation. The quantum yield is understood here as a fraction of two energies: the energy of re-emitted fluorescent photons, , over the energy of absorbed excitation photons, , where and are frequencies of fluorescent and excitation electromagnetic radiations, respectively, and are numbers of re-emitted fluorescent and absorbed excitation photons, respectively, and h is Planck’s constant.Next, we assume that fluorescent intensity propagating toward each pixel of the CCD camera results from coherently influenced waves resonantly scattered by the fluorescent particles in a voxel of the domain, which is larger than several central Fresnel zones [22, 23]. This assumption implies rectilinear propagation of the fluorescent light in accordance with Huygens’ principle. Therefore, we have
Multiply scattered fluorescent intensity, , satisfies the equation
where is the fluorescence average intensity, which is defined analogously to Eq. (2) and satisfies the Helmholtz equation: Therefore, the fluorescence intensity in the direction s is approximated by In contrast to Eq. (7), the direct fluorescent radiation is included into the expression of the observed fluorescence intensity, providing the possibility to consider fluorescent targets located near the boundary.2B. Inverse Problem: Variational Framework
The image reconstruction algorithm is based on minimization of an appropriate cost functional [24]. The cost functional is constructed from (i) the error norm (the objective function), which is the norm of difference between measured and computed intensities on the surface of an scattering object, (ii) the Lagrangian terms, and (iii) the regularization term. Minimization of the cost functional is performed in the Fourier domain due to its simplicity in contrast to the time domain minimization. Measured excitation and fluorescent intensities were transformed into the Fourier domain by applying the fast Fourier transform. The first variation of the cost functional results in a system of partial differential equations. The solution of this system is used for iterative reconstruction of optical and fluorescent parameters. Equations are solved on a Cartesian mesh by utilizing the finite volume numerical scheme [25, 26] and the ray tracing technique [27]. Although, in our experimental setup, relative positions of the laser source and the CCD camera are fixed, while the object under study is rotated, computationally we rotate source positions and the CCD camera in the opposite direction. Then, the variational problem is formulated as a minimization problem of the following cost functional:
In Eq. (14), the error norm is given as where and are experimentally recorded excitation and fluorescent intensities, respectively. The functions , , and are introduced for the sake of convenience. Thus, assuming collimated intensity measurements, we define where represents the direction of light rays coming to the CCD array, and N is the number of source–camera positions. Similarly, the functions and represent sampling in space and frequency: where M is the number of pixels of an image, L denotes the number of samples in the Fourier domain (ω), and the vector denotes the surface points visible by the CCD camera. The factor is the pixel area and, therefore, gives the total visible area. The form of is chosen in order to simplify a variational procedure. Thus, the function allows one to replace a sum over surface points visible by the CCD camera with a volume integral. Analogously, the function replaces a sum over samples in the Fourier domain with an integral.Lagrangian terms are introduced as
In Eqs. (18, 19), denotes the inner product in , J and are the adjoint excitation and fluorescent intensities, respectively, and ψ and are adjoint excitation and fluorescent average intensities, respectively.Then, we define a vector of four unknowns at every point of the scattering domain as
and introduce the regularization term In Eq. (21), are Tikhonov’s regularization parameters, and , where k denotes the iteration number starting from the initial guess .The reconstruction algorithm results from equating the first variation to zero:
Variation of functions , , ψ, and recovers Eqs. (1, 8, 11, 3, 12). Variations of I and result in equations satisfied by adjoint intensities: Note that , where n is the outward normal to the CCD array. Therefore, adjoint intensities propagate from the camera. Between an object and the camera, free space propagation is assumed. Variations of and u give where an asterisk denotes complex conjugate quantities. Finally, variations of the vector of parameters x lead to the following system: where , 2, 3, 4, and In numerical implementation of the scheme, dimensionless variables are introduced as , where are typical values of parameters, and are dimensionless variables of the order of . Because holds, a few small terms in Eqs. (28, 29) are neglected.3. Experimental Details
3A. FRET Construct Preparation and Image Acquisition
This new approach for reconstruction of time-gated tomographic fluorescence imaging data was applied to the experimental data acquired for [1]. Briefly, these data were acquired by imaging two mice with right hind legs cotransfected by electroporation with plasmids for eGFP and mCherry fluorescent proteins that were expressed separately and two mice with legs transfected with plasmids encoding an (eGFP–mCHerry) FRET construct in which eGFP (donor) was coupled to mCherry (acceptor) by a short peptide flexible linker. This transfection by electroporation primarily targets the tibialis anterior muscle, but it may include some of the surrounding muscles in the anterior lateral quadrant of the leg. For imaging, the anesthetized mice were positioned on a motorized rotation stage with the target leg held under tension approximately along the axis of rotation. The excitation beam was selected from an ultrafast supercontinuum source to have a central wavelength of and bandwidth . This was focused onto the surface of the leg and adjusted laterally such that the excitation spot (diameter ) was approximately centered on the rotation axis. The excitation and fluorescence radiation emerging from the leg was imaged onto a time-gated optical image intensifier synchronized to the excitation source and read out by a CCD camera with the intensifier’s gate width set to . The excitation and fluorescence radiation was sampled over a range at intervals and this acquisition procedure was repeated for each of 12 angular orientations of the mouse leg ( separation) and for two vertical illumination positions. For each image set, the different positions of the excitation spot were recorded to provide the source positions for subsequent tomographic image reconstruction. Following optical acquisition, the mouse was euthanized by terminal anesthetic overdose. The leg was then embedded in 1% agarose (Sigma) and imaged separately using magnetic resonance imaging (MRI). More information on the acquisition of the experimental data can be found in [1].
3B. Mesh Construction
Computational meshes of mice legs were constructed from their shadows (shadowgrams) taken by the same CCD camera, which was used for wide-field time-gating image acquisition. Shadowgrams were acquired at angular steps. They were treated as Radon transforms of completely opaque objects. Then, the filtered backprojection algorithm was applied, providing stacks of two-dimensional slices. These stacks were stored in 3D arrays and segmented with the help of the Perona–Malik method [28]. Figure 1 shows the constructed meshes. The voxel size of the resulting Cartesian meshes corresponds to a pixel size of the CCD camera. The total number of computational voxels varies from 300,000 to 500,000, depending on leg size. Although the meshes are quite dense, solution of the Helmholtz equation with consequent ray integration takes no more than a few minutes, depending on the mesh size. For instance, solution of the Helmholtz equation takes about by using the finite volume numerical discretization scheme and the conjugate gradient method, and the Siddon algorithm takes less that for mesh containing approximately 300,000 voxels on a processor.
4. Results and Discussion
Reconstruction results of four legs are presented in Fig. 2. The first and third rows [(a)–(e) and (k)–(o)] show reconstruction of optical and fluorescent parameters of two control mouse legs transfected with unlinked eGFP and mCherry. The second and fourth rows [(f)–(j) and (p)–(t)] present mouse legs transfected with FRET construct. The first column shows meshes of four mouse legs. Meshes are cut at heights where slices are taken. At these heights, the legs were expressing fluorophore and, therefore, slices contain the highest contrast in the parameter reconstructions. The second column presents reconstruction of the absorption coefficient, , while the third column displays the scattering coefficient, , in . Reconstructed quantum yield, η, and lifetime, τ, are displayed in the fourth and fifth columns. The lifetime is given in nanoseconds. The maximum values of reconstructed lifetimes are 2.97 and 2.54 for control mice legs, and 2.0 and 1.86 for legs transfected with FRET construct. Reference values obtained from cytosol preparations of lifetimes are approximately for unlinked eGFP and for eGFP linked with mCherry [1]. Therefore, reconstructions demonstrate a consistent lifetime contrast for FRETing and non-FRETing cases. Reconstructed quantum yield, however, does not present such contrast. Reconstructed values of the yield vary within a interval. Moreover, accuracy in lifetime reconstruction indicates that the error in reconstructed values can reach 15%–20% or more. Inaccuracy of reconstruction of fluorescent parameters is mostly attributed to the ill-posed nature of the inverse problem and autofluorescence present in the visible part of the light spectrum ().
Imaging in the visible part of the light spectrum presents other difficulties. The relatively high absorption coefficient precludes imaging of larger parts of an animal body. Indeed, reconstruction maps of optical parameters indicate quite high absorption of about for tissue and higher values for bones. Furthermore, strong absorption results in lower signal-to-noise ratio. Nevertheless, localization of bones (tibia) are relatively accurate. Figure 3 shows absorption and scattering coefficients merged with slices obtained by post mortem MRI without changing position by embedding the mouse leg in agar gel. Bone’s localization is slightly shifted only for the first mouse leg as a result of inaccuracy of the mesh construction procedure.
4A. Telegraph Equation Approximation
It is instructive to compare the results presented above against the TE-based reconstruction. Modifications made to the scheme in order to use it with the TE are simple ones. Thus, the light transport is described by Eqs. (3, 12) only. In the functional [Eq. (14)] the term [Eq. (18)] is dropped and the error norm is replaced with
where , u, , and are excitation and fluorescent average intensities on the surface of light scattering objects. The functions and are computed from experimentally measured intensities and , respectively. For conversion of intensities to average intensities, the method of Schwarzschild–Schuster [17] was applied in a way described earlier [19]. In the reconstruction algorithm, Eqs. (23, 26) are replaced with two Helmholtz equations for adjoint average integrated intensities, ψ and [19], and all terms containing I, , J, and are dropped in Eqs. (28, 31).Reconstruction of optical and fluorescent param eters of four mouse legs are shown in Fig. 4. The layout of the figure is the same as in Fig. 2. The last two columns, where the yield and lifetime are shown, indicate that the TE model provides very similar reconstructions of fluorescence parameters. On the other hand, reconstructions of optical parameters reveal differences. First of all, only slices (g) and (h) indicate relatively correct location of a bone. On the other slices, reconstructions of optical parameters are failed. Failure of recovering optical parameters is a direct consequence of (i) low signal-to-noise ratio present in the visible part of the light spectrum, (ii) the ill-posed nature of the inverse problem, and (iii) image mapping procedure. Thus, the conversion procedure according to the Schwarzschild–Schuster method implies an optically homogeneous medium near the boundary in addition to the assumption of smoothness of the boundary. Presence of an absorber and/or scatterer near the boundary results in inaccuracy of the image conversion procedure.
5. Summary
We described a new technology for in vivo imaging of FRET in mice that employs optical tomography and FLIM. A novel light transport model was introduced and developed. Our model combines the TE approximation with rays optics. This approximation is well adopted to intensity-based measurements of small scattering objects, which is used in wide-field time-gated data acquisition. Reconstruction results were compared against results obtained by using the TE approximation. Comparison indicates noticeable improvement in recovering of optical parameters in favor of the suggested approach, while reconstruction of fluorescent parameters shows similar results.
This work was supported by the Wellcome Trust Technology Development grant 086114.
1. J. McGinty, D. W. Stuckey, V. Y. Soloviev, R. Laine, M. Wylezinska-Arridge, D. J. Wells, S. R. Arridge, P. M. W. French, J. V. Hajnal, and A. Sardini, “In vivo fluorescence lifetime tomography of a FRET probe expressed in mouse,” Biomed. Opt. Express 2, 1907–1917 (2011). [CrossRef] [PubMed]
2. J. R. Lakowicz, in Principles of Fluorescence Spectroscopy (Plenum, 1999).
3. D. L. Andrews and D. S. Bradshaw, “Virtual photons, dipole fields and energy transfer: a quantum electrodynamical approach,” Eur. J. Phys. 25, 845–858 (2004). [CrossRef]
4. S. R. Arridge and J. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009). [CrossRef]
5. V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. 8, 1–33 (2006). [CrossRef] [PubMed]
6. M. A. O’Leary, D. A. Boas, X. D. Li, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 21, 158–160 (1996). [CrossRef] [PubMed]
7. A. T. N. Kumar, J. Skoch, B. J. Bacskai, D. A. Boas, and A. K. Dunn, “Fluorescence-lifetime-based tomography for turbid media,” Opt. Lett. 30, 3347–3349 (2005). [CrossRef]
8. A. T. Kumar, S. B. Raymond, G. Boverman, D. A. Boas, and B. J. Bacskai, “Time resolved fluorescence tomography of turbid media based on lifetime contrast,” Opt. Express 14, 12255–12270 (2006). [CrossRef] [PubMed]
9. R. E. Nothdurft, S. V. Patwardhan, W. Akers, Y. Ye, S. Achilefu, and J. P. Culver, “In vivo fluorescence lifetime tomography,” J. Biomed. Opt. 14, 024004 (2009). [CrossRef] [PubMed]
10. V. Gaind, K. J. Webb, S. Kularatne, and C. A. Bouman, “Towards in vivo imaging of intramolecular fluorescence resonance energy transfer parameters,” J. Opt. Soc. Am. A 26, 1805–1813 (2009). [CrossRef]
11. V. Gaind, S. Kularatne, P. S. Low, and K. J. Webb, “Deep-tissue imaging of intramolecular fluorescence resonance energy-transfer parameters,” Opt. Lett. 35, 1314–1316 (2010). [CrossRef] [PubMed]
12. C. Dunsby, P. M. P. Lanigan, J. McGinty, D. S. Elson, J. Requejo-Isidro, I. Munro, N. Galletly, F. McCann, B. Treanor, B. Onfelt, D. M. Davis, M. A. A. Neil, and P. M. W. French, “An electronically tunable ultrafast laser source applied to fluorescence imaging and fluorescence lifetime imaging microscopy,” J. Phys. D 37, 3296–3303 (2004). [CrossRef]
13. J. McGinty, J. Requejo-Isidro, I. Munro, C. B. Talbot, P. A. Kellett, J. D. Hares, C. Dunsby, M. A. A. Neil, and P. M. W. French, “Signal-to-noise characterization of time-gated intensifiers used for wide-field time-domain FLIM,” J. Phys. D 42, 135103 (2009). [CrossRef]
14. V. Y. Soloviev, K. B. Tahir, J. McGinty, D. S. Elson, M. A. A. Neil, P. M. W. French, and S. R. Arridge, “Fluorescence lifetime imaging by using time gated data acquisition,” Appl. Opt. 46, 7384–7391 (2007). [CrossRef] [PubMed]
15. V. Y. Soloviev, J. McGinty, K. B. Tahir, M. A. A. Neil, A. Sardini, J. V. Hajnal, S. R. Arridge, and P. M. W. French, “Fluorescence lifetime tomography of live cells expressing enhanced green fluorescent protein embedded in a scattering medium exhibiting background autofluorescence,” Opt. Lett. 32, 2034–2036 (2007). [CrossRef] [PubMed]
16. J. McGinty, V. Y. Soloviev, K. B. Tahir, R. Laine, D. W. Stuckey, J. V. Hajnal, A. Sardini, P. M. W. French, and S. R. Arridge, “Three-dimensional imaging of Förster resonance energy transfer in heterogeneous turbid media by tomographic fluorescent lifetime imaging,” Opt. Lett. 34, 2772–2774 (2009). [CrossRef] [PubMed]
17. V. V. Sobolev, A Treatise on Radiative Transfer (Van Nostrand, 1963).
18. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999). [CrossRef]
19. V. Y. Soloviev, C. D’Andrea, P. S. Mohan, G. Valentini, R. Cubeddu, and S. R. Arridge, “Fluorescence lifetime optical tomography with discontinuous Galerkin discretisation scheme,” Biomed. Opt. Express 1, 998–1013 (2010). [CrossRef]
20. V. Y. Soloviev and S. R. Arridge, “Optical tomography in weakly scattering media in the presence of highly scattering inclusions,” Biomed. Opt. Express 2, 440–451 (2011). [CrossRef] [PubMed]
21. V. Y. Soloviev and S. R. Arridge, “Fluorescence lifetime optical tomography in weakly scattering media in the presence of highly scattering inclusions,” J. Opt. Soc. Am. A 28, 1513–1523 (2011). [CrossRef]
22. H. C. van de Hulst, in Light Scattering by Small Particles (Dover, 1981).
23. M. Born and E. Wolf, Principles of Optics (Pergamon, 1968).
24. J. Nocedal and S. J. Wright, Numerical Optimization (Springer-Verlag, 1999). [CrossRef]
25. V. Y. Soloviev and L. V. Krasnosselskaia, “Dynamically adaptive mesh refinement technique for image reconstruction in optical tomography,” Appl. Opt. 45, 2828–2837 (2006). [CrossRef] [PubMed]
26. V. Y. Soloviev, “Mesh adaptation technique for Fourier-domain fluorescence lifetime imaging,” Med. Phys. 33, 4176–4183 (2006). [CrossRef] [PubMed]
27. R. L. Siddon, “Fast calculation of the exact radiological path for a three-dimensional CT array,” Med. Phys. 12, 252–255 (1985). [CrossRef] [PubMed]
28. P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Trans. Pattern Anal. Machine Intell. 12, 629–639 (1990). [CrossRef]