Modeling of Photochemical Reactions in a Focused Laser Beam.

Fluorescent materials play a prominent role in the qualitative and quantitative measurement of scientific phenomena of importance in biotechnology and biomedical applications. Photodegradation of fluorophores is a process that determines the accuracy and sensitivity of such measurements. This is the motivation for developing methods for accurately measuring fluorophore photodegradation rates. Recently, illumination consisting of short pulses has been used to examine the decay of photochemical reaction products. However, the time resolved measurements are difficult to interpret since the photodegradation process usually involves multiple time scales. The frequency domain measurement technique discussed here looks at the frequency response of a fluorescent sample to a frequency modulated illuminating light. The photodegradation rate is obtained by interpreting the frequency domain measurements in terms of traditional impedance concepts. In the measurements described in this paper, a focused laser beam is used to illuminate a sample of slowly flowing fluorescent solution. The laser beam is assumed to have a Gaussian power distribution hence illumination is spatially non-uniform in the region of interest. The photochemical reaction rates depend on power, so they will also vary with the position in the beam. However in the case of photodegradation of fluorophores, the measurement of the resulting decrease in fluorescence is given in terms of the radiation emitted from the entire illuminated region. In this work we present a mathematical description of the time evolution of the fluorescence response integrated over a non-uniformly illuminated domain. As a result of our analysis, an experimentally accessible and tractable mathematical model Eq. (19) and Eq. (30) is obtained from a more fundamental description given by Eq. (4) and Eq. (5). The model is used to create a functional form for fitting experimental measurements from a lock-in amplifier.

1. Introduction

The photodegradation of fluorophores is an important process that determines the accuracy of many biological assays [1], the efficacy of photodynamic therapy [2], and the phototoxicity of drugs. The measurement of the apparent rate of photodegradation was traditionally performed by measuring the decrease in fluorescence emission under constant illumination [3]. Recently illumination consisting of short pulses has been used to examine the decay of photochemical reaction products [4]. The time resolved measurements are difficult to interpret since the photo-degradation process usually involves multiple time scales. The frequency domain measurement technique looks at the response for each frequency of modulation of the illuminating light. Harmonic equations are simpler to solve and, thus, interpretation of frequency domain measurements can be performed in terms of traditional impedance concepts [5]. In the measurements described below, a focused laser beam is used to illuminate the sample. The power distribution of a Gaussian laser beam is strongly dependent on spatial location; therefore, the rates of photochemical reactions, which depend on power, will also depend on location in the beam. Since absorption of radiation leads to local heating and concentration gradients, convective and diffusive mass transport will be present in the vicinity of the focused beam. A slow flow is imposed on the fluorophore solution in order to dominate and thus minimize the effect due to uncontrollable convective and diffusive mass transport. In the case of photodegradation of fluorophores, measurement of the resulting decrease in fluorescence usually measures the radiation emitted from the entire illuminated region. Therefore an analysis method has to be developed which describes the time evolution of the integrated fluorescence response from a nonuniformly illuminated region. In previous work [6] a frequency domain measurement technique was interpreted with a model which assumed a uniform beam profile. This work extends the model to a laser beam with a Gaussian profile. In it we derive an experimentally accessible mathematical model Eq. (19) and Eq. (30) from first principles Eqs. (4) and (5). A dramatic simplification of these equations is a consequence of the strongly disparate time scales in the kinetics imposed by the slow rate of degradation relative to the relaxation and absorption rates The first section contains a motivation and presentation of the equations that are taken as a fundamental description of the experimental setting (also see Fig. 1). The equations are based on a two-state model with excitation and relaxation. In Sec. 2, a perturbation scheme employing the disparate time scales is used to reduce the system of two partial differential equations, describing the kinetic process to a single partial differential equation Eq. (14) for the total fluorophore population (to leading order in the perturbation parameter). The biggest reduction in mathematical complexity comes from writing the total measured fluorescence in terms of the total fluorophore population averaged over the Gaussian light beam distribution. This formula is introduced in Sec. 3, Eq. (9) and there follows the derivation of a single ordinary equation for the time evolution of the spatial average Eq. (30) that completes the model. The mathematical arguments that justify the manipulations performed in the main body of the paper can be found in the appendices. Section 4 discusses the time independent case and a functional form for fitting experimental measurements from a lock-in amplifier. The form, obtained from fitting values predicted from the model, shows good agreement with actual data. In a forthcoming paper, a harmonic approximation of the solution of the model equations will be developed. It will be possible then to express the fit parameters in terms of the physical parameters in the experiment.

Fig. 1

A graphical representation of the components included in the kinetic model discussed in the text. The rate constant k, represents photodegradation that leads to non-fluorescent species. The rate constants k and k represent absorption and relaxation to the ground state respectively. N0 and N1 represent the populations of the ground and excited states.

2. Kinetic Equations

As a starting point, consider a model for reactions at a fixed point in space. This simplified model will be used to motivate the details of the analysis. Assuming no mass transfer from surrounding fluid, Fig. 1 shows a two state model where absorption occurs from the ground state with rate k, and the excited state population of size N1, relaxes back to the ground state with a rate k. Here N0 is the population of the ground state molecules. Intact fluorophores change into non-fluorescent species with a rate k. The kinetics are summarized by a set of equations

To interpret k, Eq. (1) has to be supplemented by a model of the photochemical processes that are responsible for the conversion of the fluorophores into non-fluorescent species [6]. In the present manuscript, there will be no attempt to interpret the rate constant k. The discussion below will focus on two problems inherent in using Eq. (1) to describe real experimental conditions. The first problem is that the laser beam has a spatial intensity distribution leading to spatial variation of the absorption rate. The second problem is that convection and diffusion currents will change the local concentration of molecules in an uncontrolled manner. The second problem will be reduced significantly by introducing a flow that will change the concentration of fluorophores in a predictable way.

The absorption rate depends on the power density of the incident light. The absorption rate will be written as k = σ [7], where σ(cm2) is the molecular absorption cross section, and I is the incident photon flux (1/s cm2). During the experiment, the laser irradiance P(W/cm2) is measured. The irradiance P can be converted into a photon flux, I, by dividing P by the energy per photon. Explicitly, I = P where the conversion factor is given by . Here λ is the wavelength, n is the index of refraction, h is Planck’s constant, and c is the speed of light in vacuum. This yields the result k = σ = k. The constant k relates the laser irradiance to the molecular absorption rate.

Figure 2 shows the geometric relation of the laser beam, the cuvette holding the sample, the flow, and the detector. The laser beam is incident along the z axis which points out of the plane of the paper. The irradiance of the modulated beam that illuminates the sample will be written as:

Fig. 2

The geometric arrangement of cuvette walls (thick lines), laser beam propagating in the z direction (out of the plane of the paper), the fluid velocity in the x direction, and the detector located on the y axis. The shaded circle represents the Gaussian profile of the laser beam.

P0(x, y) is the time independent component, and ΔP(x, y) is the amplitude of the modulated component, ω is the radial modulation frequency, is the frequency in Hz, and t is time. The beam will be assumed to have a Gaussian profile [8]. where P0 is the total power (watts) of the laser beam and w is the width of the beam. The function f (x, y) is defined such that its integral over all the x, y plane is normalized to 1. The time independent component and the modulation amplitude have the same spatial dependence given by f (x, y). Clearly in this situation, the absorption rate will have a strong spatial dependence. Consequently reaction products may have concentration gradients and uneven heating may cause thermal mass transfer. The kinetic model of photodegradation, given by Eq. (1), will be extended to include a spatially dependent absorption rate and a constant flow of the solution containing the reaction product. It will be assumed that the flow dominates all mass transport (diffusive and convective) into and out of the illuminating region and provides a well defined initial concentration of fluorophore in the illuminated region.

To account for the spatial variation of absorption and the flow of fluid carrying fluorophores past the laser light, Eq. (1) is modified to,

Here k (x, y, t) = k(x, y, t). The velocity of the fluid, denoted by v is assumed to be a constant and in the direction of the x axis in Fig 2. Note in contrast to Eq. (1) the absorption rate is now a function of space and time. To completely specify the solution of these equations, the number of particles flowing into the entrance of the apparatus must be specified and the number of particles at the beginning of the observation period (t = 0) must also be given.

Equation (5) states that prior to entering the laser beam there is no change in the populations of the fluorophore states. The distance L is some length along the x axis where the laser beam intensity is effectively zero. Eq. (5) also gives the initial conditions at t = 0 when the laser beam is turned on.

3. Perturbation Analysis of Kinetic Equations

The faster processes, absorption and relaxation dominate the time variation at least initially. Returning to the model for the fixed position in Eq. (1), this suggests that after an initial transient the quasi-equilibrium condition holds,

Substitution of this condition into Eq. (1) results in a single equation for the total fluorophore population N = N0 + N1.

In this section this reduction is extended to the case of non-uniform illumination and is made rigorous. Here the coupled pair of partial differential equations, Eq. (4), will reduce to a single partial differential equation Eq. (14), accurate to first order in a perturbation expansion. Justification of the operations leading to these simplifications can be found in Appendix A.

The absorption and the excited state relaxation rates are usually of the order 108 s−1, while the value of k is less than 103 s−1. The constants in Eq. (4) can therefore be re-scaled as , where ε = 10−8 and , are of the order of 1. Equation (4) can then be rewritten as,

We will treat N0 (x, y, t,) and N1 (x, y, t) as functions of ε. Let and be the leading terms for N0 (x, y, t, ε) and N1 (x, y, t, ε) respectively in a perturbation expansion for ε > 0 and large time t.

The validity of the expansion in ε is discussed further in Appendix A. Substituting Eq. (7) into Eq. (6) and equating coefficients of equal powers of ε on both sides of the equation gives for ε = 0

This implies that as ε → 0

Thus the relation between the populations of ground and excited states holds in this approximate sense for all larger times. (This is not valid for t ≈ 0. When the illumination is turned on, there is an initial transient period where the populations are changing rapidly and therefore where the left hand side of Eq. (6) is not proportionate to ε.) Equating the coefficients of the first power of ε gives:

Letting , be the total fluorophore population in the zero order approximation, we can add the two parts of Eq. (10) and obtain

The right side of Eq. (11) depends on the zeroth order approximation to the population of excited state. Using Eq. (9) and the definition of the total fluorophore population it is possible to obtain the relation which is correct to zero order in ε. Substituting Eq. (12) into the right hand side of Eq. (11) gives the following equation for N (x, y, t) valid up to first order in ε

Equation (10) has been rewritten as a kinetic equation for the total fluorophore number valid to the first order in ε. Equation (13) can be rewritten further by setting b = k/k, η = k. Eq. (13) becomes

The boundary condition for N(x, y, t) at x = − L is inherited from those imposed on N0 and N1 at x = − L in Eq. (5)

The initial condition at t = 0 is N(x, y, 0) = N0 since no reaction takes place prior to turning on the illumination.

The last term in Eq. (14) gives the flow related flux of fluorophores into the illuminating region as a function of position. The flow of fluorophores into the volume illuminated by the laser beam is assumed to be uniform over the entire volume. The photodegradation depends on the

4. Spatial Averaging and the Fluorescent Signal

The measured fluorescence signal per unit distance in the z direction is given by where A represents the characteristics of the optical measurement instrument, krad is the rate of radiative decay of the excited state population, and N1(x, y, t) is the population of the excited state. (The radiative decay contributes to the overall relaxation rate given by k). The population of the excited state can be approximated by Eq. (12) leading to

The symbol b was defined previously. Equation (17) together with Eq. (14) constitute a description of the measurement valid to first order in ε. Equation (17) suggests that instead of finding the solution N(x, y, t) at every spatial point it may be better to find the variation of the averaged concentration defined by Eq. (18)

The average is performed with a weighting factor equal to the spatial distribution of the laser irradiance. The weighting factor provides a measure of the contribution of fluorophores at each point in space to the total fluorescence signal. It will be shown below that Eq. (14) can be converted into an equation for 〈N(t)〉 as defined by Eq. (18), and that the fluorescence signal can be written as where to a very good approximation α is a constant and P(t) = P0 + P1cos(ω t) ≡ P0 + Δ(t). For actual computations, the constant α is evaluated in terms of the time independent solution of Eq. 14. In this section, an ordinary differential equation for the time evolution of 〈N(t)〉 will be derived (Eq. (30)). The fluorescent signal can therefore be obtained (with small error) by solving Eq. (30) and then using Eq. (19). Together these equations represent an experimentally accessible mathematical model of the fluorescent signal under excitation by a focused laser beam.

To find an equation for 〈N(t)〉, first multiply Eq. (14) by 1 + bP(x, y, t) and obtain (below r represents the coordinates x,y).

Multiply Eq. (20) by f (x, y) and integrate over the beam cross section region, R. Each term will be treated in detail. To simplify the notation integrals in the x–y plane will be represented by a single integral in r = (x, y) with dr = dxdy. The left side of Eq. (20) becomes which can be expanded to

Inserting the explicit form of the power function we get

The second part of the above equation will be rewritten by introducing the quantity,

The final result for the left side of Eq. (20) after multiplication and integration is therefore

(Note that the definition given in Eq. (21) will be used in deriving Eq. (19)).

The first term on the right side of Eq. (20) after the operations becomes where is defined in Eq. (21). The second term on the right side of Eq. (20) becomes, after multiplication and integration,

To treat the second term of the above equation, we introduce a second function defined by

The second term of the right hand side of the transformed Eq. (20) now reads

In Appendix B we will show that the functions and are, up to a small error, independent of t and can be replaced by constants α and α′ respectively. The error introduced into equations Eq. (21), Eq. (23), and Eq. (25) by this approximation is discussed in Appendix C.

Collecting terms, Eq. (20) can be written after multiplication by f (x, y) and integration, up to small varying functions of t, as where

The final step in the derivation of an equation for 〈N(t)〉 consists of finding an approximation to the weighted average spatial derivative in Eq. (27). For small widths w the major contribution to the integral comes from a region in the vicinity of the beam which we define by 0 ≤|x|,|y|≤ κ w where we may choose the constant κ, 0 ≤ κ ≤ 3. If we replace the partial derivative in the integrand by a backwards difference quotient with step h = − κ w the value of N at the left endpoint of the interval is N0 (to a good approximation). Thus for w small enough we have, where o(1)is a term that tends to zero as w → 0. To guarantee an error that is smaller than the dominant terms of Eq. (26) we have to require that w is sufficiently small and κ is large. In fact experimental conditions govern the actual choice of w. The following heuristic argument governed the choice of κ.

We observed that in the vicinity of the origin (x = 0, y = 0) the spatial derivative of the time independent solution, N(x, y) discussed below, has a functional form which is almost identical to the beam distribution function f (x, y) as shown by the two overlapping curves in Fig. 4. The difference between − f (x, y) and ∂N(x, y)/∂x (both normalized by dividing by the maximum of their respective absolute values) is magnified 50 times and shown in Fig. 4 by the upper curve. The weighted average of the un-normalized beam distribution function π · width2 f (x, y) is equal to 0.5. Therefore the weighted average of the spatial derivative will be approximately equal to 0.5 times the value at the origin

Fig. 4

The correspondence between the Gaussian beam profile function − f (x, y), shown by the solid line, and the steady state function , shown by the dotted line. Both functions are normalized by their maximum values in order to emphasize functional similarities. The difference between the two functions, magnified 50 times, is shown by the curve with two valleys. The graphs in Fig. 4a are for y value of 0 cm, and the graphs in Fig. 4b are for y value of 0.005 cm.

The derivative at the origin can be approximated in the usual way by the difference of the function N(x, y, t) evaluated at two points spanning the origin and divided by the distance between the two points. However for the purpose of this analysis a more useful approximation uses the difference 〈N (t)〉 − N0 divided by an effective distance which reproduces the value of the derivative. For the case of the Gaussian beam, the best estimate of the derivative at the origin was found (using the time independent solution) to be (〈N (t)〉 − N 0)/w*0.637. This yields the following approximate expression for the average spatial derivative given in Eq. (28).

Using Eq. (29), the approximate kinetic equation for 〈N (t)〉 becomes where Δ(t) = P1 cos(ωt).

5. Validation of the Model Equations and Their Derivation

This section discusses the time independent form of Eq. (30) (see Eq. (31). Since an analytical form of the solution exists for the time independent case of Eq. (14), direct comparisons of the fluorescent signal obtained from the time independent version of Eq. (30) are possible. In both cases the fluorescence signal is obtained from Eq. (19) (using the steady state solutions of Eq. (14) and Eq. (30). A second important reason for considering the time independent solution is that the spatial average in Eq. (18) is dominated by the spatial average of the time independent solution itself. As a consequence, the time varying functions α (t) and α′ (t) can be replaced by constants calculated from the time independent solution (see B 9 in Appendix B).

The time independent form of Eq. (14) corresponds to a condition of constant illumination, constant flow, and constant initial fluorophore concentration.

Figure 2 shows the geometrical arrangement of the flow channel, the laser beam, and the detector. The laser beam is along the z direction, the detector along the y axis, and the flow takes place along the x axis. The channel dimensions are d y = 0.4 cm along the y axis, d z = 1 cm along the z axis, and 5 cm along the x axis. The velocity of the fluid in the center of the channel is set to 0.025 cm/s. No uncertainties are given since these values will be used to model the response of the system. The width of the Gaussian laser beam is set to w = 50 · 10−4 cm. The laser power is P0 = 0.03 W, the vacuum wavelength is λ = 488 · 10−7 cm. The number of photons per Joule, Pc, is calculated by Pc = λ n/hc = 3.314 · 1021 J−1, where n = 1.35 is the index of refraction of water and hc is the product of Planck constant and the speed of light. The extinction coefficient is set to ε = 82000 L/cm mol, yielding a cross section σ = 1000ε/0.4343 N = 3.135 · 10−16 cm2 where N is the Avogadro’s number. The average absorption rate is given by , and the relaxation rate is given by the inverse of the measured lifetime of the excited state, k = 1/τ = 3.33 · 108 s−1. The photodegradation rate k is set to 170 s−1. Using these values it is possible to calculate the parameters that enter into Eq. 31,b = σ P/k = 3.117 · 10−6 cm2 s/J and η = k = 5.3 · 10−4 cm2/J.

The solution of Eq. 31 is written as where L = 0.02. Equations (21) and (24) can be used to calculate α = 6324 and α′ = 8462, respectively, using Eq. (32). Note that A is a constant which is set to 1. Figure 3 shows the resulting N(x, y) for the parameters given above. The solution given by Eq. (32) can be used to predict the time independent fluorescence signal by substituting it into Eq. (19). The result can be compared to fluorescence calculated from the steady state of the approximate kinetic Eq. (30), i.e.,

Fig. 3

(a) The steady state solution of Eq. (31) given by Eq. (32). The physical parameters that enter into the equation are discussed in the text. The function N(x, y) is normalized to 1. The laser beam is incident along the normal to the x-y plane. As expected, the number of intact fluorophores decreases with increasing laser beam intensity. The decreased concentration of intact fluorophores persists along the x axis due to the imposed flow. (b) The steady state x-derivative of the function N(x, y). The derivative is finite along the profile of the laser beam.

The two values of the fluorescent signal (FS) agree to within 1 % up to 0.1 W incident power. The approximate fluorescence signal given by Eq. (30) and Eq. (33) is almost indistinguishable from the result given by Eq. (32). It is expected that the time dependent solutions will also be described correctly by the approximate Eq. (30) which together with Eq. (19) will be taken as the fundamental description of the photodegradation process occurring in the focused laser beam. Finally the steady state solution was used to find the effective distance to use in the estimate of the average spatial derivative given by Eq. (29). The exact value of the average spatial derivative and the estimated value using Eq. (29) are − 3.932 cm− 4 and − 3.906 cm−4 respectively. Thus Eq. (29) should provide a reasonable estimate of the average spatial derivative even in the time dependent case.

The formal time dependent solution of Eq. (30) can be written as (see for example http://tutorial.math.lamar.edu/terms.aspx) where

For the case of low power, the approximations bα P(t) « 1, bα′ P(t) « 1, and a laser beam power given by P(t) = P0 + P1 cos (ω t) yield an approximate solution given by Eq. (36).

Figure 5 shows the time dependent solution at three values of the laser power for a modulation frequency of Hz. The other parameters were taken to be the same as in the case of the steady state solution. The salient features of the time dependent solution are a transient at t = 0 and a modulation of the intact fluorophore concentration at the driving frequency of the laser. As expected the average level of intact fluorophores decreases with increasing average power, and the amplitude of the modulation increases with increasing modulation amplitude. The approximate time dependent solution given by Eq. (36) was used to obtain the fluorescence signal given by Eq. (19). The fluorescence signal was detected by a computer model of a lock-in amplifier consisting of an in-phase and quadrature multipliers followed by low pass filters. The low pass filters were implemented by Fourier analyzing the signal, setting to zero all coefficients corresponding to non zero frequencies, and then performing an inverse Fourier transform. The solid circles in Fig. 6a show the ratio of the detected quadrature and in-phase components at different modulation frequencies. Most of the parameters used in Eq. (36) were the same as in the steady state calculation except for width = 10 µm, P0 = 50 mW, and k = 600 s−1. The solid line in Fig. 6a shows a fit to the calculated points of a function given by R(f) = (a* f)/(b + f2), where f is the modulation frequency and a and b are fit parameters. Clearly the calculated response and R(f) form an excellent match. The solid points in Fig. 6b show experimentally measured values of the ratio of quadrature and in-phase response. The solid line in Fig 6b is a fit to the function R(f). The match is reasonable. Thus there is a strong indication that the measured response can be represented by the function R(f), and that the fit parameters provide information about the various physical constants of the experiment. To elucidate the dependence of the fit parameters on the physical constants, it is best to use analytical solutions to Eq. (30) and obtain the analytical form of the response function R(f). This task is carried out in a forthcoming paper.

Fig. 5

The time dependent numerical solutions of Eq. (30). Each curve corresponds to a different power of the laser beam. The initial value is set to 1 in all cases. After an initial transient, the average concentration of intact fluorophores oscillates with the modulation frequency. The amplitude of the oscillation increases with increasing laser power.

Fig. 6

(a) The time dependent solutions, such as shown in Fig 5, were inserted into Eq. (19) to model the detected fluorescence signal. The fluorescence signal was inserted into a digital model of a two phase lock-in amplifier (discussed in the text). The ratio of the outputs of the quadrature and in-phase lock-in amplifiers is shown by the solid circles. The solid line is a best fit to a function R(f) = (a* f)/(b + f2), where f is the modulation frequency and a and b are parameters. (b) The solid circles show the measured response of a flowing fluorescein solution in a focused laser beam with a wavelength of 488 nm. The solid line gives the best fit to the function R(f) given in (a).

We close this section with a brief description of the mathematical justification for the time dependent case. The reduction of the original system Eq. (4) and Eq. (5) to a single partial differential equation depended on the use of a multiple time scale perturbation expansion. The validity of the technique used for this problem (matched asymptotic expansions) is discussed in Appendix A. By spatially averaging the partial differential equation we derive an ordinary differential equation for the averaged reactant concentration as well as a simplified expression for the total fluorescence. These simplifications are possible because the varying functions (Eq. (21) and Eq. (24) respectively) are nearly constant. The proof of this fact can be found in Appendix B. Finally the equation for the spatial average, as in Eq. (30), involves these functions and it is shown that replacing them by constants creates a small error. This is shown in Appendix C.

6. Conclusion

The dynamics of the populations of excited and unexcited fluorophores are described by a pair of first order partial differential equations that incorporate the kinetics of the fluorophore reaction driven by a periodically varying laser beam and the convective fluid flow in a measurement apparatus. We obtain a dramatic reduction in the mathematical complexity of the model by taking advantage of the disparate time scales for excitation, relaxation and degradation. This is followed by averaging the reactant concentration over the laser beam, which is assumed to have a Gaussian profile. A single ordinary differential equation, Eq. (30), is obtained describing the evolution of the spatially averaged reactant concentration. The measured fluorescent signal, Eq. (19), is expressed in terms of the spatially averaged reactant concentrations. These equations provide a prediction of the response whenever the photochemical reaction degrades the fluorophores flowing past the laser beam.

Mathematical justification for the approximations used in the course of the derivation is given and some comparisons to the time independent case are described. The predicted frequency response (based on the time varying solution of Eq. (30)) is used to derive a functional form for the ratio of the in-phase and quadrature parts of the fluorescent signal. This same form provides a good fit for experimental data. The fit and the connection with physical parameters of the experiment will be discussed in a forthcoming paper.