WO1997008538A1 - Fluorescence lifetime-based imaging and spectroscopy in tissues and other random media - Google Patents

Abstract

A system and method non-invasive biomedical optical imaging and spectroscopy with low-level light is described. The technique consists of a modulated light source (120) coupled to tissue (100) of a patient to introduce excitation light. Fluorescent light emitted in response to the excitation light is detected with sensor (148). The AC intensity and phase of the excitation and detected fluorescent light is provided to a processor (160) operatively coupled to sensor (148). Processor (160) employs the measured re-emission kinetics of excitation and fluorescent light to 'map' the spatial variation of one or more fluorescence characteristics of the tissue (100). The fluorescence characteristic may be provided by exogenous contract agents, endogenous fluorophores, or both. The variations is determined by solving frequency domain diffusion equations at a number of designated points in the tissue as part of a recursive estimation algorithm. Processor (160) generates an imaging signal in accordance with the spatial variation of the fluorescence characteristic for provision to an output device (164). The output device (164) displays an image corresponding spatial variation of the fluorescence characteristic which corresponds to tissue (100) to aid in the detection and diagnosis of disease.

Description

FLUORESCENCE LIFETIME-BASED IMAGING AND SPECTROSCOPY IN TISSUES AND OTHER RANDOM MEDIA

BACKGROUND OF THE INVENTION

The present invention relates to spectroscopic imaging of heterogeneous light scattering media, and more particularly, but not exclusively, relates to in vivo imaging of biologic tissue by mapping a fluorescence characteristic ofthe tissue through the detection of light emitted in response to excitation light from a time- varying light source.

The early detection of disease promises a greater efficacy for therapeutic intervention. In recent years, non-invasive techniques have been developed which have improved the ability to provide a reliable and early diagnosis of various afflictions by detecting biochemical changes in the tissue of a patient. For example, Magnetic

Resonance Imaging (MRI) has successfully monitored the relaxation of spin states of paramagnetic nuclei in order to provide biomedical imaging and biochemical spectroscopy of tissues. Unfortunately, the complexity and expense of MRI diagnostics limit its availability — especially as a means of routine monitoring for disease. Another powerful analytical technique with an increasing number of applications in the biological sciences is fluorescence spectroscopy. These applications include biomedical diagnostics, genetic sequencing, and flow cytometry. To date, there are several industrial and academic institutions developing fluorescent and phosphorescent compounds for observing pertinent metabolites and environmental conditions, such as Ca^"*^"4", pH, glucose, pθ2, and pCO2- With the development of dyes and photodynamic fluorescent agents which excite and re-emit in the near-infrared red (NIR) wavelength regime, non-invasive detection of diseased tissues located deep within tissues may also be possible since red excitation and re-emission light can travel significant distances to and from the tissue-air interface (See Wilson et al.. Time- Dependent Optical Spectroscopy and Imaging for Biomedical Applications. 80 Proceedings IEEE pp. 918-30 (1992)). As exemplified by U.S. Patent Nos. 5,421,337 to Richards-Kortum et al. and

5,452,723 to Wu et al., several investigators have suggested various procedures to differentiate diseased and normal tissues based on fluorescence emissions through non¬ invasive external measurements or minimally invasive endoscopic measuring techniques. Unfortunately, these procedures generally fail to provide a viable spatial imaging procedure. One reason imaging based on fluorescence has remained elusive is that meaningful relational measurements of fluorescence characteristics from a random, multiply scattering media, such as tissue, are difficult to obtain. For example, fluorescent intensity, which is a function ofthe fluorescent compound (or fluorophore) concentration or "uptake," is one possible candidate for imaging; however, when this property is used in an optically dense medium, such as a particulate (cell) suspension, powder, or tissue, the local scattering and absorption properties confound measured fluorescent intensities.

Besides intensity, other properties of selected fluorophores such as fluorescent quantum efficiency and lifetime are also sensitive to the local biochemical environment. As used herein, "fluorescent quantum efficiency" means the fractional number of fluorescent photons re-emitted for each excitation photon absorbed or the fraction of decay events which result in emission of a fluorescent photon. "Fluorescent lifetime," as used herein, is defined as the mean survival time ofthe activated fluorophore or the mean time between the absorption of an excitation photon and re-emission of a fluorescent photon. Like intensity, measurement of these fluorescence characteristics is often limited to well-defined in vitro applications in the research laboratory or in flow cytometry where issues such as scattering, absorption, and changing fluorophore concentrations can be controlled or measured. Moreover, these limitations generally preclude meaningful fluorescence-based imaging of hidden tissue heterogeneities, such as tumors or other diseased tissue regions which cannot be detected by visual inspection.

Thus, a need remains for a technique to non-invasively image multiply scattering tissue based on one or more fluorescence characteristics which does not require extensive information about intrinsic optical properties ofthe tissue, and takes advantage of the contrast capability offered by fluorescence yield and lifetime characteristics to aide in the identification of tissue heterogeneities. The present invention satisfies this need.

SUMMARY OF THE INVENTION

The invention relates to spectroscopic imaging of heterogeneous, light scattering materials. Several aspects ofthe invention are novel, non-obvious, and provide various advantages. While the actual nature ofthe invention covered herein can only be determined with reference to the claims appended hereto, certain features which are characteristic ofthe present invention are described briefly as follows.

One feature of the present invention is a technique for imaging a heterogeneous light scattering material. This process includes exposing the surface of a material to light from a light source and detecting an emission in response. A spatial variation of a fluorescence characteristic ofthe material is determined as a function ofthe emission with a processor. The spatial variation may be characterized by a set of values representative of the fluorescence characteristic as a function of position. An image is generated in accordance with the spatial variation which corresponds to the heterogeneous composition ofthe material. This technique may be applied in vivo to biologic tissue using extemal or endoscopic instrumentation to detect heterogeneities indicative of disease. The technique may include the introduction of a fluorescent contrast agent into the material. The fluorescence characteristic detected may be fluorescence lifetime, fluorescence quantum efficiency, a fluorophore absoφtion coefficient, fluorescent yield (a function of fluorescent quantum efficiency and fluorophore absoφtion), or another fluorescence characteristic known to those skilled in the art.

In another feature ofthe present invention, the spatial variation of a light scattering material with a heterogeneous composition is determined by establishing an estimate ofthe optical property or fluorescence characteristic variation, determining a calculated emission from the material as a function ofthe estimate, and comparing the calculated emission to a detected emission to determine a corresponding error. The estimation ofthe variation is modified, the calculated emission re-determined with this modified estimate, and the comparison repeated until the error reaches a desired minimum. An image ofthe material is generated from the modified estimate which corresponds to the heterogeneous composition.

Accordingly, it is one object ofthe present invention to map a fluorescent property of a light scattering material which varies with the heterogeneous composition ofthe material to generate a corresponding image.

It is another object ofthe present invention to provide a spectroscopic technique for non-invasively monitoring fluorescent properties of hidden tissue volumes in a living organism and to monitor selected metabolites of an organism in vivo. Yet another object is to provide a fluorescence imaging system and process to identify diseased tissue using endogenous or exogenous fluorophores as a contrast agent. This contrast may be provided by a concentration, lifetime, or quantum efficiency difference between normal and diseased tissues.

It is still another object ofthe present invention to provide an imaging technique and algorithm which is based on contrast of an optical property that is independent of local fluorophore concentration.

Further objects, features, aspects, and advantages ofthe present invention will become apparent from the drawings and description contained herein.

BRIEF DESCRIPΗON OF THE DRAWINGS

FIG. 1 is a schematic illustration of a system of one embodiment ofthe present invention. FIG. 2 is a flow chart of a process utilizing the system of FIG. 1.

FIG. 3 is a schematic representation of a tissue phantom arrangement used to demonstrate various aspects ofthe present invention.

FIGS. 4-7 graphically depict selected properties of equations used in the present invention. FIGS. 8 and 9 graphically depict convergence of simulated determinations of the spatial variation of fluorescent yield and lifetime, respectively, utilizing one embodiment ofthe present invention.

FIGS. 10-14 are images obtained from experimental examples 1-3 ofthe present invention. FIG. 15 is a schematic illustration of a system of an alternative embodiment of the present invention.

DESCRIPΗON OF PREFERRED EMBODIMENTS

For the puφoses of promoting an understanding ofthe principles ofthe invention, reference will now be made to the embodiment illustrated in the drawings and specific language will be used to describe the same. It will nevertheless be understood that no limitation ofthe scope ofthe invention is thereby intended. Any alterations and further modifications in the described device, and any further applications ofthe principles ofthe invention as described herein are contemplated as would normally occur to one skilled in the art to which the invention relates. FIG. 1 depicts system 110 of the present invention for fluorescent imaging of tissue 100. Tissue 100 has surface 101 and a heterogeneous composition as represented by regions 102, 103 underlying surface 101. Heterogeneities 102, 103 are generally not detectable by visual inspection of surface 101.

System 110 includes modulated light source 120 to supply an intensity modulated excitation light of predetermined frequency and wavelength to tissue 100 via optic fiber 123. Preferably, source 120 is a laser diode of conventional design with a modulated output in the 1-500 MHz frequency range and a monochromatic output in the 100 to 1000 nm wavelength range. The specific wavelength is selected to excite a designated fluorophore in tissue 100. Beam splitter 126 may be employed to direct a small portion ofthe excitation signal to reference sensor 128 for processing puφoses. System 110 also includes detection subsystem 140 which has optic fibers 143 to detect photons emitted from tissue 100 from a number of corresponding detection sites. Subsystem 140 includes one or more emission sensors 148. Detection subsystem 140 also includes an interference filter to obtain a selected emission wavelength corresponding to emission of a designated fluorophore in tissue 100. In one embodiment, subsystem 140 includes a single sensor 148 and the signals from fibers 143 are multiplexed. Preferably, sensors 128, 148 are Photo-multiplier Tubes (PMTs) or photodiodes but other sensor varieties, such as image intensifiers and charge-coupled devices, are also contemplated. Sensors 128, 148 and source 120 are operatively coupled to heterodyne subsystem 130. Subsystem 130 is configured to obtain information about the phase, AC, and DC intensity of light detected with sensor 128 relative to light detected with the sensor 148 using conventional laser heterodyning techniques. In one embodiment, heterodyne subsystem 130 includes a signal synthesizer phase-locked to the repetition rate of a laser used for source 120. For this embodiment, subsystem 130 includes an amplifier to gain modulate sensors 128, 148 at a harmonic of a laser repetition rate (when a pulsed laser is used) or at the modulation frequency (when a modulated laser diode is used) plus an offset to provide the desired heterodyning. In one variation of this embodiment, an 80 MHz pulsed laser repetition rate is divided down to 10 MHz and input to the synthesizer, and a heterodyning offset of 100 kHz is input to the amplifiers for sensors 128, 148.

Sensors 128, 148 are operatively coupled to processor 160. Processor 160 includes input/control device 162, output device 164, and memory 166. Processor 160 may be an electronic circuit comprised of one or more components. Similarly, processor 160 may be comprised of digital circuitry, analog circuitry, or both. Also, processor 160 may be programmable, an integrated state machine, or a hybrid combination thereof. Preferably, input device 162 is a keyboard or input control of a conventional variety, and output device 166 is a Cathode Ray Tube (CRT) based video display, printer, or other image display system known to those skilled in the art. Memory 166 is preferably ofthe electronic (e.g. solid state), magnetic, or optical variety ofthe type readily available for use with electronic controllers or processors. Furthermore, Memory 166 may include an optical disk memory (CD), electromagnetic hard or floppy disk media, or a combination of these. FIG. 2 depicts one mode of operation of system 110 as process 210. Process

210 includes mapping the spatial variation of fluorescence yield and lifetime with processor 160 and generating an image signal in accordance with the map. Output device 164 is configured to display an image in response to the image signal. Process 210 begins by introducing a fluorescent contrast agent into tissue 100 in stage 212. This agent provides a source of fluorescent emission for detection by subsystem 240. The configuration ofthe modulated light source 120, heterodyne subsystem 130, and detection subsystem 140 is designed to accommodate the excitation and emission properties ofthe selected fluorescent agent. In other embodiments, endogenous fluorophores may be alternatively or additionally employed and system 110 adjusted accordingly.

In stage 214, tissue 100 is excited by light source 120 configured according to the selected fluorophore. In stage 216, the phase, θ_ODS, and log of AC intensity, M_0DS, ofthe emission at each detection site "i" relative to the excitation light from source 120 are determined at the heterodyne (or offset) frequency. For "Di" number of detection sites, the detected or observed phase and AC intensity are indexed by "i" using the following notation: (θ₀b_s)i and (M_0DS)i, respectively. Processor 160 stores the relative phase and AC intensity information in memory 166.

In stage 218, a two dimensional grid is established for an area of tissue 100 selected for imaging, and a matrix of grid points is established and indexed by "j". A uniform seed value for the fluorescent yield, y^; - ("n ax-^rn j' ^an(^ *^^e fl^uorescent lifetime, (τ)j, at each grid point j is assigned. These values are an initial homogeneous guess of the yield and lifetime values which are modified in later stages. The term "η" is the quantum efficiency ofthe fluorophore which varies with the environment ofthe surrounding ofthe fluorophore. The term " Hax-»m" ^s *^^{e aDSO}Φ^uon coefficient for the fluorophore and is the product ofthe extinction coefficient ofthe fluorophore based on the natural log and the concentration ofthe fluorophore. As a result, the yield, y = ημ_aχ→m, is influenced by the surrounding metabolism and the uptake ofthe fluorophore. The uptake of certain known fluorophores vary with the type and condition of host tissue, providing another fluorescence characteristic useful to detect disease. The contrast provided by these properties is largely independent of fluorophore concentration. The initial estimate of fluorescent yield and lifetime are stored in memory 166 by processor 160 for later use. After establishing this initial estimate ofthe fluorescence characteristics of yield, ημa_χ_»rn' ^anc^ ^^e^^mei^ ^τ> processing loop 220 is entered in stage 230. Preferably, the stages of processing loop 220 are executed by processor 160 via preprogrammed software, dedicated hardware, or a combination of both as appropriate. To aid in understanding various mathematical aspects of process 210 and loop 220, the following table of selected variables is listed:

c velocity of light;

D(r) optical diffusion coefficient; Di number of detection sites; f modulation frequency;

I identity matrix; i detection site index;

J Jacobian matrix relating the sensitivity at each grid point, j, to the response at each detection site; j grid point index;

J_j?i individual elements ofthe Jacobian matrix J; k source index;

M log of AC intensity of modulated fluorescent light position; m index to multiple modulation frequencies; n average index of refraction; r position (in two or three dimensions);

Sk number of modulated light sources;

S(r, ω) source term for the modulated light at position r and frequency ω; ^'

Greek χ² merit function representing the least squares error; Φ_x(r, ω) complex number representing photon flux in the frequency domain at position r and frequency ω; η quantum efficiency of fluorescent probe or dye; μ_a average absoφtion coefficient; absoφtion coefficient ofthe fluorescence light by both the non- fluorescing chromophores and fluorophore;

^a_x absoφtion coefficient ofthe excitation light by both the non- fluorescing chromophores and fluorophore;

^μaχ→c adsoφtion coefficient due to non-fluorescing chromophores;

^ x→.m adsoφtion coefficient of excitation light by fluorophores; μ's effective scatting coefficient; θ phase-shift of one modulated light wave to another; τ lifetime of activated probe or dye at location r; ω angular modulation frequency, given by 2πf;

Subscripts obs observed or experimental data; x excitation light; and m fluorescence or emission light.

In stage 230, phase and relative AC intensity at each detection site "i" is calculated as a function ofthe initial estimates of yield and lifetime for each grid point j. The calculated phase and intensity are represented at each detection site i as (θ_m)i and (M_m)j, respectively. The values for (θ_m)j and (M_m)j are determined using the diffusion equation approximation ofthe radiative transport equation. The diffusion equation approximation describes the spatial and temporal transport of light in tissues or multiply scattering media. A coupled frequency domain diffusion equation can be used to predict the excitation and emission fluence rates, Φ_x(r, ω) and Φ_m (r, ω), respectively, at any location r within the selected grid of tissue 100 via equations (1) and (2):

V • [D_x(r)VΦ_x(r, ω)]-[μ_aχ(r)+i ω/cJΦxCr, ω) +S_x(r,ω) = 0 (1)

V • [D_m(r)VΦ_m(r, ω)]-[μ_am(r)+i ω/c_n]Φ_m(r, ω) +S_m(r,ω) = 0 (2)

The source term for the excitation light S_x(r,ω) is due to the sinusoidally modulated light at an angular frequency ω = 2ωf where f is typically in the MHz frequency range. The first term in both ofthe diffusion equations (1) and (2) represents the diffusive or "random-walk" transport of light where D_{x m} is the optical diffusion coefficient of equation (3) as follows:

and μ_a and μ'_s are the absoφtion and isotropic scattering coefficients, respectively, for tissue 100, the medium of interest. The optical properties are dependent on the wavelength of light and thus are different for the excitation light from source 120 (subscript x) and fluorescent emission detected with subsystem 140 (subscript m). The total absoφtion coefficient at the excitation wavelength, μ_aχ, is due to contributions from non-fluorescing chromophores as well as from fluorophores responsive to the excitation wavelength. The total absoφtion coefficient is given by the sum of absoφtion coefficients due to non-fluorescing chromophores, μ_aχ→c, and fluorophores ax-_Mrr Generally it may be assumed that the absoφtion experienced at the fluorescent wavelength is due primarily to non-fluorescing chromophores. The velocity of light in tissue is (^c/n where n is the average index of refraction. The source term for the fluorescent emission is dependent on the excitation light fluence, Φ_x(r, ω) and is given by equation (4) as follows: S_m(r,ω) = ημa_χ→m (r)Φ_x(r, ω)[(l - iωτ(r))/(l + crA(r)²)] (4)

This term arises from the Fourier transform ofthe fluorescence decay term in the time domain following an incident pulse of excitation light where: τ is the fluorophore lifetime, η is the quantum efficiency, and the absoφtion coefficient, a_x_»_m> *^s ^ product ofthe extinction coefficient based on natural log and the concentration ofthe fluorophore in the ground state. As previously indicated, the combined product, ημ_aχ→._m, is termed the fluorescent yield, y, and is proportional to the generated fluorescence fluence. Substitution of equation (4) into equation (2) facilitates determination of Φ_m for each grid point "j." The solution ofthe diffusion equations (1) and (2) for the two-dimensional area defined by the grid points "j" may be readily extended to three dimensions to estimate spatial variation of one or more fluorescence characteristics in a selected volume with "r" corresponding to position in three dimensions.

Both diffusion equations (1) and (2) are linear complex elliptic equations that can be solved as boundary value problems for the complex quantities Φ_χ(r, co) and Φ_m(r, ω). This solution employs the method of finite differences to create corresponding finite difference equations. These difference equations are utilized to obtain an approximate solution at each grid point, j. This method of solution is described in other contexts in Fulton et al., Multigrid Method for Elliptic Problems. A Review, 114 American Meteorological Society pp. 943-59 (May 1986); and B.W. Pogue et al., Initial Assessment of a Simple Svstem for Frequency Domain Diffuse Optical Tomography. 40 Physics in Medicine and Biology pp.1709-1729 (1995). One preferred method of performing this solution is with the MUDPACK routines described in Adams, J.C, MUDPACK: Multigrid Portable Fortran Software for the Efficient Solution of Linear Elliptic Partial Differential Equations. 34 App. Math Comp. p.133 (1989). For the solution ofthe diffusion equations, it is assumed that Φ_m,χ(r,co) = 0 on the surface 101 of tissue 100 which is known as the zero fluence boundary condition. It should be recognized that other boundary conditions may be selected and the method of solution varied accordingly.

The diffusion equations (1) and (2) may be solved for a complex number for Φ_m at each grid point, j. The detected signal at the surface is proportional to the normal component ofthe gradient ofthe photon fluence. To approximate the signal at detector site "i" located on surface 101 of tissue 100, the Φ_m value at an internal grid point closest to the site is selected which follows from the relationship that the normal component ofthe photon fluence gradient is proportional to Φ_m just inside surface 101. The calculated phase-lag, θ_m, and the log of AC intensity, M_m, at the detection sites "Di" are calculated from the imaginary and real parts ofthe complex Φ_m value with respect to the phase and the AC intensity of source 120.

The diffusion equations (1) and (2) provide insight into the sensitivity of changing the fluorescent optical properties of tissue 100 on θ_m and M_m measured at the detector sites i. This insight results from a series of calculations while fixing various parameters ofthe diffusion equations (1) and (2). These calculations assume circular tissue phantom 300 with an embedded, heterogeneity 302 hidden in phantom background 303 as illustrated in FIG. 3. A two-dimensional grid is established for phantom 300 and may easily be expanded to three dimensions. Under these simulated conditions, a large value is assigned to absoφtion coefficients for both excitation and fluorescent light at all grid points outside the simulated tissue phantom. The four sources S1-S4 of FIG. 3 (Sk = 4) are simulated by assigning an arbitrary complex number at a grid point near the surface closest to each source. The twenty detection sites D1-D20 of FIG. 3 (Di = 20) are simulated by using the calculated values determined from Φ_m at the grid point "j" closest to the detection site. The simulated solutions to diffusion equations (1) and (2) were obtained in two dimensions for a 65 x 65 grid covering a 100mm diameter circular tissue phantom 300 with a circular, embedded heterogeneity having a 30mm diameter and located at the center ofthe tissue phantom 300 (this location differs slightly from the configuration of heterogeneity 302 of FIG. 3). The simulated measurements of fluorescent phase-shift and AC intensity are reported for 20, equally spaced, circumferentially located detection sites D1-D20. The modulation frequency, f, was set equal to 150 MHz. The optical properties ofthe heterogeneity and the background are shown in Table 1 as follows:

Table 1

In order to evaluate the influence of ημ_aχ→m> θm a d M_m were computed at each detection site D1-D20 as the value of ημ_aχ→._m i the heterogeneity increased from 10^-4 mm^-1 to 10^_1 mm^-1 and as ημ_aχ→m i the background 303 was maintained constant. The lifetime, τ, was set equal to 1 ns for both the object and the background causing contrast due to differences in ημ_aχ→m- The plots of θ_m and M_m are shown in FIGS. 4 and 5 respectively for one active source SI. As ημ_aχ→._m of heterogeneity 102 increases to higher values, the AC intensity approaches an upper limit similar to what is expected in dilute non-scattering solutions. FIG. 5 shows how the fluorescent phase- shift, θ_m, decreases as the absoφtion coefficient due to the fluorophore, μ_a "•_vΛ"¹ .ill is decreased 10 to 100 times the background. From these simulations, M_m appears to be directly dependent upon changes in ημa_x→m of a simulated tissue heterogeneity 102 whereas θ_m is indirectly dependent on ημ_aχ→m due to changes in photon migration. In order to evaluate the influence of τ, θ_m and M_m were calculated at each detection site D 1 -D20 as the values of τ in the heterogeneity varied from 10" * ns to 10³ ns and the value of τ in the background was held at 1 ns. The background ημ_a was set to lO^'^rnn and η a_x→m f°^r *h^e heterogeneity was set to 10"3mm"l. As shown in Fig. 6, the detected AC intensity increases as τ decreases. Fig. 7 illustrates the values ofthe fluorescent phase-shift at each detection site as the lifetime ofthe heterogeneity is changed from 0.1 ns to 1000 ns. At a given modulation frequency (150 MHz in this calculation), θ_m first increases, reaches a maximum and then subsequently decreases as τ is increased from 0.1 ns to 1000 ns. Therefore, both θ_m and M_m at each detection site D1-D20 appear to be directly influenced by the value of lifetime in the heterogeneity.

Referring back to FIG. 2, in stage 240, the calculated emission phase and intensity, (θ_m)j and (M_m)j, are compared to the measured emission phase and intensity, (θ₀bs)i and (M_0DS)i, for each detection site "i" to identify a difference or "error" between the measured and calculated values. Because (ημ_aχ→._m)j impacts (M_m)j , this comparison is posed in the form of the merit function χμ2 of equation (5) as follows:

Sk Di χ_μ2 = (l/Sk) Σ (1/Di) Σ [((M₀bs)i - (M_m)i) σM]2 (5) k=l i=l

where σ^ is the typical standard deviation of noise in M_m, taken to be 0.01; Sk = number of excitation source sites indexed by k; and Di = number of detection sites indexed by i. The goal ofthe algorithm is to minimize χ„2 by appropriate updates of iμax— » _m)j- After ^an initial update of (ημ_{a →} )j another merit function in terms of (τ): participates in the comparison of stage 240. This merit function, χ_τ2 , is presented as equation (6) as follows:

Sk Di χ_τ2 = (l/Sk) Σ (1/Di ) Σ [((Mo_bs)i - (M_m)i)/σM] + [((θ_obs)i - (θ_m)i)/σ_θ]2 (6) k=l i=l

where σø is the typical standard deviation of noise in (θ_m)i, taken to be 1 degree; Sk = number of excitation source sites indexed to k; and Di = number of detection sites indexed to i. Since the lifetime influences both (θ_m)i and (M_m)j, the phase and AC intensity are used in equation (6). After the comparison of stage 240 is performed by calculating the merit functions χμ2, χ_χ2 ₁ control flows to conditional 250 to test whether the comparison of

the observed values, (θ₀b_s) and (Mob_s)i, to the calculated values (θ_m)j and (M_m)_j

via the merit functions meets a selected convergence criteria. This criteria corresponds

to the degree of tolerable error in determining the yield and lifetime values. For one

embodiment, convergence is achieved when any ofthe following three quantities, (i)

χ2, (ii) change in χ2 in successive iterations of loop 220, (ii) relative change in χ2 in

successive iterations of loop 220 is less than a predetermined threshold value of 1.0 x 10"2. In other embodiments a different comparison calculation and associated

conditional may be employed as would occur to one skilled in the art. If conditional

250 is satisfied, control flows to stage 270 and loop 220 is exited; however, if the

criteria is not satisfied, execution of loop 220 continues in stage 260.

In stage 260, the yield, (y)^; = Cημa_χ-»_m i" ^an(i life ime, (τ)j, for each grid point j is updated so that these values may reach the minimum error corresponding to the

comparison stage 240 and conditional 250 test. In order to update these values, Jacobian matrices are used which describe the sensitivity ofthe response at each

detection position i to changes in (y)j = (η aχ->m)j» ^anc* life*"¹¹⁶' (^τ)j> ^{at eacn} 8" point, j. Three Jacobian matrices are employed: J (M, ημ_aχ__{> m}); J(M, τ); and J(θ,

τ). The elements Ji j of these Jacobian matricies are given by Ji j =

[θM_{/(5(ημ_aχ→m)j]; Jij = [δMj/dτj]; and Ji = [θθ/θτj], respectively. These elements may be calculated by solving the diffusion (1) and (2) four times for each grid point, j

to obtain M_m i and θ_m,j calculated with (τ )j and (τ+ δτ)j and with (ημa_x-_→._m)j an (ημ_a + δημ_aχ→m)j. From least squares rninimization, the update to yield and lifetime is calculated. In one preferred embodiment, this updating algorithm is adapted from an algorithm used to reconstruct images obtained by electrical impedance tomography like the algorithm suggested by Yorkey, et al, Comparing reconstruction

Algorithms for Electrical Impedance Tomography. 34 Transactions in Biomedical

Engineering pp. 843-52 (1987). The Jacobian Matrices are used to solve for update

vectors, [Δημ_aχ→._m ] and [Δτ], to estimated yield and lifetime vectors, [ημ_{aχ→ m} ] and [τ], respectively. These vectors are of a dimension corresponding to the number of

grid points. At each iteration through loop 220, the following Jacobian equations (7)

and (8) are solved to determine the update for the estimated yield and lifetime vectors:

M_JJ, and M_m are the observed and calculated vectors of the log of AC intensity at each ofthe i detection sites, respectively, θ_m and θ_m are the observed and calculated vectors ofthe phase lag at each ofthe i detection sites, respectively. Due to the ill-conditioned nature ofthe Jacobian matrices, the terms

are added as part of a Marquardt rrmώnization scheme where I is an identity matrix. The parameters λ_j or λ2 are adjusted via a Maquardt-Levenberg type algorithm ofthe type disclosed in Press et al., Numerical Recipes: The Art of Scientific Computing. (Cambridge University Press, 1992). Conventional numerical methods are employed to solve the simultaneous linear algebraic equations resulting from the Jacobian matrix equations (7) and (8). The Jacobian matrices are re-calculated at each iteration through loop 220. It has been found that equations (7) and (8) provide a way to select appropriate changes to the yield and lifetime estimates; however, other numerical approaches to recursively iterate to acceptable estimates as would occur to one skilled in the art are also contemplated. Once the update is complete, control returns to stage 230.

If the convergence criteria is satisfied in conditional 250, then estimation of yield and lifetime for the grid points has reached an acceptable minimum and control flows to stage 270. In stage 270 an image signal is generated by processor 160 from the spatial variation ofthe yield and/or lifetime fluorescence characteristics. This image signal is sent to output device 164 which displays an image in response. Because the fluorescence characteristics of yield and lifetime typically vary with the biologic environment ofthe fluorophore, this image is generally indicative of tissue variation and offers the capability to detect heterogeneities 102, 103.. For example, laser diodes capable of supplying Near infrared (NIR) light that can penetrate tissue several centimeters, and fluorescent contrast agents responsive to NIR light may be used to provide a viable imaging system. In one embodiment, this system is adapted for use with an endoscope. Besides yield and lifetime, the spatial variation of other fluorescence characteristics useful to distinguish diseased tissues may be mapped using the diffusion equations (1) and (2). Such alternative fluorescence characteristics include, but are not limited to, quantum efficiency η and/or fluorescent absoφtion coefficient μa_x_>_m determined as separate properties independent ofthe yield product. In another embodiment ofthe present invention, the photon fluence equation and Jacobian estimation process is adapted to determine a map of a designated fluorophore uptake concentration. For this embodiment, a first map of chromophore adsoφtion coefficients μa_x→c and scattering coefficients μ'_s are determined in the absence ofthe designated fluorophore by estimating the chromophore adsoφtion coefficient μ_aχ→c and scattering coefficient μ'_s at each grid point j in place ofthe yield and lifetime estimates. Diffusion equation (1) for Φ_x(r, ω) may be used in conjunction with modified Jacobian equations (7) and (8) to create this first map. The modification substitutes the chromophore adsoφtion and scattering coefficients in place ofthe yield and after adaptation to accommodate these new characteristics as follows:

and

where ri_j = Sk and n^ = Di. After generating the first map, the designated fluorescent contrast agent is introduced, and the total adsoφtion coefficient μ_aχ is determined by substituting μ_aχ in place of μ_aχ→._c equations (9)-(l 1) to obtain a second map ofthe total adsoφtion coefficient. Noting that μ_aχ = μ_aχ→._m ⁺ μa_χ_ c' ^an^ ^{t ιat} ^^e UP*²^ °f the fluorescing contrast agent is directly proportional to μa_x→m' up*⁸^ concentration may be mapped by deteπnining a difference between the adsoφtion coefficient variations for the first and second maps. This "difference map" may then be used to generate an image corresponding to the uptake concentration.

Another alternative embodiment measures the emission responsive to each of a number of light source modulation frequencies f. The total number of different frequencies employed is designated Mf. To obtain this additional data, an iteration of loop 220 is performed for each frequency f indexed to m. The number of sources, Sk and detection sites Di are indexed to k and i, respectively. This additional data may be used to enhance imaging results obtained with system 110 or to permit reduction ofthe number of detection sites or excitation source sites in the evaluation. A representative merit function corresponding to this additional data is given in equation (12) as follows:

χ_τ ² =

Mf Sk Di

(1/Mf) Σ (1/Sk) Σ (1/Di) Σ [((M_obs)i - (M_m)i)/σ_M]2 + [((θ_obs)i - (θ_m)i)/σ_θ]² (12) m=l k=l i=l

Besides fluorescence yield and lifetime, the multi-frequency method can be employed to map other optical characteristics of interest. Besides a sinusoidally modulated light source, the present invention may be adapted to operate with a pulsed or other time- varying excitation light source in alternative embodiments.

FIG. 15 depicts an optical system 410 of another embodiment ofthe present invention. This system includes modulated light source 420 with laser driver 422, operatively coupled laser diode 424, and reference frequency generator 426. Source 420 is configured to deliver modulated light to tissue phantom 400, and the re-emitted light from the phantom is focused onto a gain modulated image intensifier 430 via 50mm lens 432. Intensifier 430 includes a photocathode face which converts photons to electrons, a Multi-Channel Plate (MCP) which multiplies the electronic signal by avalanche multiplication, and a phosphorescent screen which converts electrons into an optical image. Preferably, intensifier 430 is a fast intensifier, ofthe variety manufactured by Litton Electronics, Inc., which enables modulation by applying a DC bias and an RF signal from amplifier 428 between the photocathode and the MCP. For this example, the modulation ofthe image from intensifier 430 is phase-locked to the laser diode 424 by a 10 MHz output signal from synthesizer 426. By modulating the laser diode 424 and the image intensifier 430 at the same frequency, a steady-state image results on the phosphor screen. U.S. Patent No. 5,213,105 to Gratton et al. provides additional background concerning certain aspects of this technique. The image from the phosphor screen is focused through interference filter 433 on a Charge Coupled Device (CCD) camera 434 via 150 mm macro lens 436. Camera 434 has a 512x512 array of CCD detectors configured to provide a corresponding pixelated image. Camera 434 is operatively coupled to processor 460 of a similar configuration to processor 160 previously described.

Following each acquired image, a phase delay between the image intensifier 430 and the laser diode 424 is induced by stepping the phase ofthe image intensifier 430 to values between 0 and 360 degrees with the frequency synthesizer 452 under the control of processor 460. Since the gain modulation of image intensifier 430 and laser diode 424 occurs at the same frequency, homodyning results in a steady phosphorescent image on intensifier 430 which is dependent upon phase. Preferably, control between synthesizer 452 and processor 460 is obtained by a conventional GPIB interface.

Images from the phosphorescent screen ofthe image intensifier 430 are then gathered at each phase delay. The incremental phase delayed images are then used to generate a map of phase-shift and intensity modulation ratio between the excitation and re-emitted light from phantom 400. By applying interference or appropriate optical filters, the emission light may be selectively separated from the excitation light and measured. Camera 434 output may be processed by processor 460 using process 210. The present invention will be further described with reference to the following specific examples 1-3. It will be imderstood that these examples are illustrative and not restrictive in nature. Examples 1-3 involve the computer simulation of the process 210. Simulations of this kind, including the simulation of tissue, are an acceptable means of demonstrating fluorescent spectroscopic imaging performance to those skilled in the art. The examples use simulated values obtained by solving the diffusion equations (1) and (2 ) for θ_m and M_m under the conditions of table 2 as follows:

Case ^a„_c ^„ μ_Sχ or τ ^"π _a Gaussian Gaussian

(background) (background) Noise Noise in log of AC in phase

(mm"¹) (mm^-1) (mm"¹) (ns) (mm^*1) intensity σø (degrees)

5.1 0.0 0.0 1.0 10.0 1.0 x IO-⁵ 0.01 0.1

5.2 1.0 x 10-3 0.0 1.0 10.0 1.0 x IO"⁵ 0.01 0.1

5.3 0.0 0.0 1.0 10.0 1.0 x 10-5 0.01 1.0

Table 2

The examples simulate tissue phantom 300 of FIG. 3 having a 100mm diameter. Values of θ_m and M_m were computed at each of the D1-D20 detection sites of FIG. 3 in response to the 4 modulated light sources S1-S4 located at the periphery. The excitation light modulation frequency f was simulated at 150 MHz. Diffusion equations (1) and (2) were solved to provide 80 simulated values of θ_m and M_m correspondmg to the various combinations of detection and source sites (Sk * Di = 4x20 = 80). Gaussian noise with a standard deviation of 0.1 degrees (or a liberal 1 degree) in θ_m and 1% in M_m were superimposed on the diffusion equation solutions. Adapted MUDPACK routines were used to solve the diffusion equations (1) and (2) on a SunSparclO computer. These obtained data sets were used as simulated input data to process 210 for examples 1-3. The results are shown in tables 3 and 4 are as follows:

Case Area, object 1 Location, object 1 Area, object 2 Location, object 2

(mm²) (x,y), (mm, mm) (mm²) (x,y), (mm, mm)

5.1 706.0 (expected) (60,60) (expected) not applicable not apphcable

742.2 (obtained) (60.8,58.5) (obtained)

5.2 706.0 (expected) (60,60) (expected) not apphcable not applicable

703.1 (obtained) (59.4,58.3) (obtained)

5.3 314.1 (expected) (32.3,67.7) (expected) 314.1 (expected) (67.7,32.3) (expected) i 381.0 (obtained) (34.0,67.7) (obtained) 342.0 (obtained) (65.0,35.0) (obtained)

Table 3

Table 4 Example 1

Example 1 reconstructs fluorescent yield and lifetime with no absoφtion due to non-fluorescing chromophores. To simulate the experimental data for this example, the fluorescent yield, iμa_χ-»._m)j_> f°^r ^ background and the heterogeneity 302 were chosen as 1 x 10"^ mm"' and 1 x 10~³ mm"* respectively and the fluorescence lifetime, (τ );, for the background and the heterogeneity 302 chosen as 10 ns and 1 ns respectively. During the execution of loop 220, no a priori knowledge of either the heterogeneity 302 location or the background fluorescence properties was assumed and a uniform guess of lx 10^"^ mm^-1 and 10 ns was given for the fluorescence yield, iμax- »._m)j' ^an<^ lifetime, (τ );, respectively. Convergence was achieved in less than 50 iterations of Loop 220 (computational time on a SunSparclO: 2 hours) for a two dimensional 17 x 17 grid. The average values of ημ_aχ_v_m and τ in the grid points which occupy the simulated object converge within 50 iterations to ημ_aχ→m ⁼ 0- 3 x 10"³ mm^-1 and τ =1.03 ns are illustrated in FIGS. 8 and 9, respectively. FIGS. 10 and 11 illustrate the reconstructed images from the mapped values of ημ_aχ_ _m and τ, respectively, and are representative ofthe expected images. The images were smoothed by inteφolation in examples 1-3 to remove spurious points which had unphysically high values, but were surround by values within a physically achievable range. These spurious values were replaced by the average background fluorescence yield and lifetime obtained from simulation of loop 220.

The average values of ημ_aχ→. _m in the grid points which occupy the simulated background converge within 50 iterations to 9 x 10"^ mm^" . The value ofthe background converges to 5.4 ns. The dependence ofthe final images on the choice of the initial guess was examined by providing an initial uniform guess of 1 x 10^"4 mm^"1 and 10ns for 0iμa_χ->_m)j' nd lifetime, (τ)j, respectively. This resulted in similar images to those obtained in FIGS. 10 and 11.

The location of heterogeneity 302 was identified as consisting of all the grid points with ημ_aχ→. _m higher than 35% (arbitrarily chosen) ofthe peak value ofthe ax— »m (FIG. 10). The average of the coordinates of all the identified object grid points was the position (60.8, 58.5) which is close to position (60, 60) that was used to simulate the experimental data. As listed in Table 3, the area ofthe heterogeneity based upon our arbitrary definition for identification was 72 mm2, close to that used to generate our simulated experimental data.

Example 2

Example 2 reconstructs fluorescent yield and lifetime with a simulated chromophore absoφtion configured to mimic tissue. The same hidden heterogeneity as well as optical parameters and simulation equipment were used as described in Example 1 except that a uniform background chromophore absoφtion coefficient, μ_aχ→ of 1 xlO"³ mm"¹ was used to generate the simulated experimental data. While excitation light propagation was not employed for image reconstruction, we considered this optical property known to estimate the best possible performance for inverse image reconstruction under physiological conditions. The two-dimensional reconstructed spatial map ofthe fluorescence yield, (ημ_{aχ→ m})j, and lifetime, (τ);, are shown in FIGS. 12 and 13, respectively. As shown in Table 3, the mean value of location ofthe object according to our criterion based on ημax→m occurred as position (59.4, 58.3) consistent with the conditions used to simulate the experimental data. The dimension ofthe heterogeneity based upon our arbitrary definition for identification (all grid points with Wax-*™ higher than 35% ofthe maximum) were 703 mm^ which is close to that used to generate our simulated experimental data. The average values of ημ_a and τ in the grid points which occupy the simulated object converge within 50 iterations to the values of ημ_aχ→._m ⁼ °-⁸ xlO^"3 mm^'1 and τ = 0.7 ns consistent with the values used to generate the simulated experimental data (see Table 3). The average values of ημ_{a →} and τ in the grid points which occupy the simulated background converge within 50 iterations to values similar to that reported for Example

1. Example 3 Example 3 simulated two hidden heterogeneities in the tissue phantom (not shown in FIG. 3). In this case, the same optical parameters were used as described in example 1 except that the fluorescence yield ημa_x_i._m f°^r the objects 1 and 2 was chosen as 1 x 10"³mm^-1 and 2 x lO'^mm"¹ respectively and lifetime τ for the heterogeneities chosen as 1 ns and 2 ns, respectively. A 33 x33 grid was employed instead of a 17 xl7 grid. An image corresponding to the mapping of yield is depicted in FIG. 14.

All publications and patent applications cited in this specification are herein incoφorated by reference as if each individual publication or patent application were specifically and individually indicated to be incoφorated by reference. While the invention has been illustrated and described in detail in the drawings and foregoing description, the same is to be considered as illustrative and not restrictive in character, it being understood that only the preferred embodiment has been shown and described and that all changes and modifications that come within the spirit ofthe invention are desired to be protected.

Claims

CLAIMSWhat is claimed is:

1. A method of imaging, comprising the steps of: (a) exposing a surface of a light scattering material to an excitation light from a light source, the material having a heterogeneous composition under the surface;

(b) detecting a fluorescent emission from the material in response to step (a);

(c) establishing an estimate of spatial variation of a fluorescence characteristic ofthe material; (d) determining a calculated emission as a function ofthe estimate;

(e) comparing the calculated emission to the emission detected in step (b) to determine an error;

(f) providing a modified estimate of spatial variation ofthe fluorescence characteristic and repeating steps (d) through (f) until the error reaches a desired minimum; and (g) generating an image ofthe material from the modified estimate, the image corresponding to the heterogeneous composition ofthe material.

2. The method of claim 1, further comprising introducing a fluorescent contrast agent into the material.

3. The method of claim 1, wherein step (f) includes applying a Jacobian matrix.

4. The method of claim 1 , wherein the expected emission is determined as a function of a diffusion equation.

5. The method of claim 1, wherein step (e) includes comparing intensity and phase ofthe calculated emission to intensity and phase ofthe emission detected in step (b).

6. The method of claim 1, wherein the material includes a fluorescent contrast agent and the fluorescence characteristic is a function of at least one of fluorescence quantum efficiency, fluorescence lifetime, and concentration ofthe fluorescent contrast agent.

7. The method of claim 6, wherein the fluorescence characteristic is a function of fluorescence quantum efficiency, the light source is intensity modulated at a predetermined frequency, step (e) includes comparing AC intensity and phase ofthe calculated emission to intensity and phase ofthe emission detected in step (b), step (f) includes applying a Jacobian operator, and the expected emission is determined as a function of a photon fluence rate within the tissue.

8. A method of imaging, comprising:

(a) exposing a living biologic tissue to an excitation light from a light source;

(b) detecting a fluorescent emission from the tissue in response to said exposing; (c) determining a spatial variation of a fluorescence characteristic within the tissue as a function ofthe emission with a processor; and

(d) generating an image ofthe tissue in accordance with the spatial variation.

9. The method of claim 8, further comprising introducing a fluorophore into the tissue.

10. The method of claim 8, wherein the fluorescence characteristic is a function of at least one of fluorescence lifetime, fluorescence quantum efficiency, or fluorescence absoφtion.

11. The method of claim 9, wherein said exposing includes positioning the source adjacent a surface ofthe tissue, and said detecting includes sensing emitted light at the surface.

12. The method of claim 8, wherein said exposing includes positioning a number of modulated light sources adjacent a surface ofthe tissue and said detecting includes sensing emitted light at a number of locations along the surface.

13. The method of claim 8, wherein said determining includes (i) establishing an estimate of spatial variation ofthe fluorescence characteristic, (ii) determining a calculated emission as a function ofthe estimate, (iii) comparing the calculated emission to the fluorescent emission of said detecting to determine an error, (iv) providing a modified estimate of spatial variation ofthe fluorescence characteristic, and (v) repeating said determining and said comparing with the modified estimate until the error reaches a desired minimum.

14. A method of imaging, comprising: (a) introducing a fluorescent agent into a light scattering biologic tissue with a heterogeneous composition;

(b) exposing a surface ofthe tissue to light from a light source to excite the agent;

(c) detecting an emission of light from the tissue in response to said exposing;

(d) determining a number of values representative of a fluorescence characteristic of the tissue as a function of position and the emission; and

(e) generating an image in accordance with the values, the image corresponding to the heterogeneous composition ofthe tissue.

15. The method of claim 14, wherein the fluorescence characteristic is a function of fluorescence quantum efficiency.

16. The method of claim 14, wherein the fluorescence characteristic corresponds to fluorescence lifetime.

17. The method of claim 14, wherein said detecting includes sensing fluorescent light emitted in response to said exposing at a number of locations along the surface.

18. The method of claim 14, wherein said exposing includes providing light from a number of sources spaced along the surface.

19. The method of claim 14, wherein said exposing includes exciting the fluorescent agent in the tissue with modulated light at a number of different frequencies, and the values are determined as a function ofthe different frequencies.

20. The method of claim 14, wherein said determining includes (i) establishing an estimate ofthe values, (ii) determining a calculated emission as a function ofthe estimate, (iii) comparing the calculated emission to the emission of said detecting to determine an error, (iv) providing a modified estimate ofthe values, and (v) repeating said determining and said comparing with the modified estimate until the eπor reaches a desired minimum.

21. A system for imaging a light scattering tissue having a heterogeneous composition and containing a fluorophore, comprising:

(a) a light source adapted to excite the fluorophore; (b) a sensor configured to provide a detected light signal corresponding to a fluorescent emission from the tissue in response to light from said source;

(c) a processor operatively coupled to said sensor and responsive to said detected light signal to provide a number of values representative of a fluorescence characteristic ofthe tissue as a function of position, the fluorescence characteristic coπesponding to at least one of fluorescence lifetime, fluorescence quantum efficiency, and fluorescence absoφtion, said processor being configured to generate an image signal as a function of said values; and

(d) an output device responsive to said image signal to provide an image coπesponding to the heterogeneous composition ofthe tissue.

22. The system of claim 21 , further comprising a number of modulated light sources.

23. The system of claim 21 wherein said fluorescence characteristic corresponds to fluorescent yield.

24. The system of claim 21, wherein said sensor is configured to detect said emission at a number of locations along a surface ofthe tissue.

25. The system of claim 21, wherein said processor is configured to determine said values

from a comparison of a calculated emission to an observed emission derived from said detected light signal, said calculated emission being determined as a function of an estimated

spatial variation of said fluorescence characteristic, said estimated variation being updated

and said comparison being repeated until a difference between said calculated emission and

said observed emission reaches a desired minimum.

26. The system of claim 21, wherein said source includes a laser diode and said sensor includes a CCD camera.