Förster Resonance Energy Transfer imaging in vivo with approximated Radiative Transfer Equation

We describe a new light transport model that we have applied to 3-D image reconstruction of in vivo fluorescence lifetime tomography data applied to read out Förster Resonance Energy Transfer in mice. The model is an approximation to the Radiative Transfer Equation and combines light diffusion and rays 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.


Introduction
We have 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 fluorescence lifetime imaging (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 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 biomolecualr interactions in live mice by introducing a new technique for tomographic image reconstruction of weakly scattering samples that utilises a novel, computationally inexpensive, approximation to the Radiative Transfer Equation (RTE) that replaces the TE.The performance of these two light transport models are compared for experimental data acquired from fluorescence protiens expressed in mouse legs [1].Briefly, the hind leg anterior tibialis muscle of live mice were transfected by electroporation and in vivo time-resolved images of the donor (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 succesful.
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, compuationally 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 the 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 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 (CT), 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.
The paper is organized as follows.Mathematical and computational details are described in the next section, where the first subsection is devoted to the direct problem of light transport, while the second subsection describes the image reconstruction approach.Experimental details and mesh construction procedure are outlined in Section 3. In the last section we compare reconstruction results based on the suggested model with results obtained by using the TE.

Europe PMC Funders Author Manuscripts
Europe PMC Funders Author Manuscripts

Methodology 2.A. Direct problem
Light transport in scattering media is modelled 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 telegraph equation approximation and takes into account the probability of the quantum exit of a photon from the medium in a given direction directly [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: (1) 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, μ s , and the absorption coefficient, μ a .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 telegraph equation (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 whole solid angle leads to the Helmholtz equation, which is the Fourier image of the TE: ( where the Helmholtz differential operator is given by: the diffusion coefficient κ is expressed in terms of the attenuation coefficient as: (4) In the source term u 0 represents the average intensity of the direct excitation radiation I 0 (s), which satisfies the equation (5) Europe PMC Funders Author Manuscripts Europe PMC Funders Author Manuscripts where q 0 is the amplitude and r 0 is the location of a laser source, which injects excitation photons into the medium in the direction s 0 .Thus, the direct excitation radiation is given by: (6) 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, r 0 , and the observation point, r, are connected by a line r = r 0 +sl.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: (7) where l max 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 I 0 (r, s 0 ) 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: (8) where (9) In Eqs. ( 9) η denotes the quantum yield and τ is the lifetime, which 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, hν f N f , over the energy of absorbed excitation photons, hν e N e , where ν f and ν e are frequencies of fluorescent and excitation electromagnetic radiations, respectively, N f and N e are numbers of re-emitted fluorescent and absorbed excitation photons, correspondingly, and h is the Planck 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: Europe PMC Funders Author Manuscripts Europe PMC Funders Author Manuscripts (10) Multiply scattered fluorescent intensity, I′, satisfies the equation: (11) where u′ is the fluorescence average intensity, which is defined analogously to Eq. ( 2) and satisfies the Helmholtz equation: (12) Therefore, the fluorescence intensity in the direction s is approximated by: (13) 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.

2.B. 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 L 2 -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 (FFT).
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: (14) In Eq. ( 14) the error norm is given as: (15) where I E and I F are experimentally recorded excitation and fluorescent intensities, respectively.The functions ς(ω), χ(r) and ξ(s) are introduced for the sake of convenience.
Thus, assuming collimated intensity measurements, we define: Europe PMC Funders Author Manuscripts Europe PMC Funders Author Manuscripts (16) where s n represents the direction of light rays comming to the CCD array, and N is the number of source-camera positions.Similarly, the functions χ(r) and ς(ω) represent sampling in space and frequency: (17) where M is the number of pixels of an image; L denotes the number of samples in the Fourier domain (ω); the vector r m denotes the surface points visible by the CCD camera.
The factor a p 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 χ(r) 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 L 2 ; J and J′ are the adjoint excitation and fluorescent intensities, respectively; and ψ and ψ′ are adjoint excitation and fluorescent average intensities, correspondingly.
Then, we define a vector of four unknowns at every point of the scattering domain as: (20) and introduce the regularization term: (21) In Eq. ( 21) α j are Tikhonov's regularization parameters, and , where k denotes the iteration number starting from the initial guess k = 0.
The reconstruction algorithm results from equating the first variation to zero: (22) Variation of functions J E , J F , ψ, ψ′ recovers Eqs.(1,8,11,3,12).Variations of I and I′ result in equations satisfied by adjoint intensities: (23) Europe PMC Funders Author Manuscripts Europe PMC Funders Author Manuscripts (24) Note that −s n = n, 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 u′ and u give: (25) (26) where asterisk denotes complex conjugate quantities.Finally, variations of the vector of parameters × lead to the following system: (27) where j = 1, 2, 3, 4, and

3.A. FRET construct preparation and image acquisition
This new approach to reconstruction of time-gated tomogrphic fluorescence imaging data was applied to the experimental data acquired for reference [1].Briefly, this data was acquired by imaging two mice with right hind legs co-transfected by electroporation with plasmids for enhanced Green Fluorescence Protein (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 (TA) muscle but may include some of the surrounding muscles in the anterior lateral quandrant of the leg.For imaging the anaesthetized mice were positioned on a motorised rotation stage with the target leg held under tension approximately Europe PMC Funders Author Manuscripts Europe PMC Funders Author Manuscripts along the axis of rotation.The excitation beam was selected from an ultrafst supercontinuum source to have a central wavelength of 480 nm and bandwidth 40 nm.This was was focused onto the surface of the leg and adjusted laterally such that the excitation spot (diameter ~ 0.5 mm) 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 synchronised to the excitation source and read out by a CCD camera with intensifier's gatewidth set to 1 ns.The excitation and fluorescence radiation was sampled over a 10 ns range at 250 ps intervals and this acquisition procedure was repeated for each of 12 angular orientations of the mouse leg (~ 30° 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 MRI.More information on the acquisition of the experimental data can be found in ref. [1].

3.B. 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 1° angular steps.They were treated as Radon transforms of completely opaque objects.Then, the filtered backprojection algorithm was applied providing stacks of 2D slices.These stacks were stored in 3D arrays and segmented with help of the Perona-Malik method [28].Figure 1 show 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 that a few minutes depending on the mesh size.For instance, solution of the Helmholtz equation takes about 75 seconds by using the Finite Volume numerical discretization scheme and the Conjugate Gradient method and the Siddon algorithm takes less that 3 seconds for mesh containing approximately 300,000 voxels on a 2.16 GHz processor.

Results and Discussions
Reconstruction results of four legs are presented in Figure 2. 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 highest contrast in the parameter reconstructions.The second column presents reconstruction of the absorption coefficient, μ a , while the third column displays the scattering coefficient, μ s , in mm −1 .Reconstructed quantum yield, η, and lifetime, τ, are displayed in the fourth and fifth columns.The lifetime is given in nanoseconds (ns).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 cytocol preparations of lifetimes are approximately 2.61 ns for unlinked eGFP and 2.13 ns 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 (0.6, 0. 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 0.1 mm −1 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 Magnetic Resonance Imaging (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.

4.A. The 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 E , u, u F and u′ are excitation and fluorescent average intensities on the surface of light scattering objects.The functions u E and u F are computed from experimentally measured intensities I E and I F , 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.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 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.

Summary
We described a new technology for in vivo imaging of Förster Resonance Energy Transfer in mice that employs optical tomography and FLIM.A novel light transport model was introduced and developed.Our model combines the Telegraph Equation 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 Telegraph Equation approximation.Comparison indicates noticeable improvement in recovering of optical Europe PMC Funders Author Manuscripts Europe PMC Funders Author Manuscripts   The lifetime is given in nanoseconds.Reconstruction results obtained with the TE approximation.
numerical implementation of the scheme dimensionless variables are introduced as , where are typical values of parameters, and are dimensionless variables of the order of O(1).Because holds, a few small terms in Eqs.(28-29) are neglected.
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.Second and fourth rows ((f)-(j) and (p)-(t)) 8) 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 Europe PMC Funders Author Manuscripts Europe PMC Funders Author Manuscripts nature of the inverse problem and autofluorescence present in the visible part of the light spectrum (480 nm).
(23)(24)(25)(26) are replaced with two Helmholtz equations for adjoint average integrated intensities, ψ and ψ′,[19] and all terms containing I, I′, J and J′ are dropped in Eqs.(28)(29)(30)(31).Reconstruction of optical and fluorescent parameters of four mouse legs are shown in Figure 4.The layout of the Figure is the same as in Figure 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.

Figure 1 .
Figure 1.Meshes of four mouse legs used in experiments.Legs (a) and (c) were transfected with unlinked eGFP and mCherry fluorophores (control mice).Legs (b) and (d) were transfected with FRET-construct, in which eGFP and mCherry are coupled by a flexible linker.

Figure 2 .
Figure 2.Reconstruction results using the approximated radiative transfer model.The first column display meshes of four legs with cuts at heights where slices are taken.The first and third rows show reconstruction of parameters of legs transfected with unlinked eGFP and mCherry.The second and third rows display parameters of legs transfected with FRET construct.The second and third columns show the absorption and scattering coefficients in mm −1 .The fourth and fifth columns display the quantum yield and lifetime reconstructions.

Figure 3 .
Figure 3. Reconstructed absorption and scattering coefficients overlayed with MRI slices.The first row (a-d) shows absorption coefficients for four mouse legs.The second row (e-h) shows scattering coefficients.Red indicates peak values.
Soloviev et al.Page 12 Appl Opt.Author manuscript; available in PMC 2012 November 08.
Europe PMC Funders Author ManuscriptsEurope PMC Funders Author Manuscripts