Research Papers: General

Fluorescence in two-photon-excited diffusible samples exposed to photobleaching: a simulation-based study

[+] Author Affiliations
Imre B. Juhász

Pázmány Péter Catholic University, Faculty of Information Technology and Bionics, Práter utca 50/A, Budapest 1083, Hungary

Árpád I. Csurgay

Pázmány Péter Catholic University, Faculty of Information Technology and Bionics, Práter utca 50/A, Budapest 1083, Hungary

J. Biomed. Opt. 20(1), 015001 (Jan 09, 2015). doi:10.1117/1.JBO.20.1.015001
History: Received August 18, 2014; Accepted December 2, 2014
Text Size: A A A

Open Access Open Access

Abstract.  We created a simulation model to investigate the characteristics of fluorescence in two-photon-excited samples. In the model, the sample is a diffusible solution of fluorophore molecules, which is divided into cubic cells and illuminated by a train of focused laser pulses described as a Gaussian beam. Simulating the state transitions according to a multilevel photodynamic model (also including photobleaching and intersystem crossing), the simulator provides the expected number and the spatial distribution of emitted photons over time. Our simulations demonstrated how the illumination laser power, diffusion, and the photodynamic parameters of the fluorophore affect fluorescence. We revealed the unusual fluorescent profile that evolves as photobleaching progresses: the most photons are not emitted from the focus (where a “dark hole” appears) but from an ellipsoid around the focus. The model could be adapted to several fluorescent techniques (such as two-photon microscopy and fluorescence recovery after photobleaching). Furthermore, it might help to optimize the operating parameters of the measurement devices (e.g., in order to reach higher image quality and lower photobleaching).

Fluorescence is a two-phase process during which photon absorption is followed by photon emission. Its high sensitivity and specificity1 are harnessed in sensing, imaging, and manipulation. In fluorescent imaging and measurement techniques, the fluorophore molecules of the sample are excited by illumination, then the fluorescent photons emitted by the relaxing molecules are collected and registered. The number of fluorescent photons refers to the number of fluorophore molecules within the excited region. Thus, based on the recorded data, a computer can depict the fluorophore distribution as a microscopic image of the sample. Moreover, the temporal changes and fluctuations of the fluorescent signal carry information about the diffusion of fluorophore molecules into and from the excited volume, which is exploited in fluorescence correlation spectroscopy (FCS)2,3 and fluorescence recovery after photobleaching (FRAP) (also known as fluorescence photobleaching recovery) measurements.4 There are different special techniques to increase the spatial resolution of imaging. In two-photon microscopes, focused light pulses generated by a mode-locked laser excite the fluorophore molecules of the sample.5 At the focal region, the intensity of illumination is high enough to provoke two-photon excitation, during which two photons are absorbed simultaneously while the molecule attains an excited state. However, moving further away from the focus, the probability of two-photon absorption decreases rapidly, because its rate is proportional to the square of the irradiating light intensity (see Sec. 2.2). The fluorescent photons can, therefore, be considered to originate from a small (subfemtoliter6) volume around the focus. The recently worked out super-resolution microscopy techniques circumvent the diffraction limit of the light microscopes formulated by Abbe and Rayleigh. These methods apply for instance special illumination patterns (as in the case of structured illumination microscopy) or special nonlinear excitation (in stimulation emission depletion microscopy), or limit the number of simultaneously emitting molecules to promote their individual localization (in photoactivated localization microscopy, fluorescence photoactivation localization microscopy, and stochastic optical reconstruction microscopy).7,8

In the illuminated sample, not only fluorescence (i.e., “excitation followed by photon emission”) takes place, since the fluorophore molecules can also undergo further possible state transitions: the excited molecule can return to the ground state via internal conversion without photon emission, it can reach the triplet state by intersystem crossing, or it can absorb more photons. These absorbed photons can cause photochemical reactions called photobleaching, which destroys the fluorescent ability of the fluorophore molecule irreversibly. Photobleaching has a significant influence on fluorescence: on one hand, it is a major limit of image quality in two-photon microscopy,9 and it also causes considerable difficulties in FCS;2,3 on the other hand, it is a fundamental process in FRAP measurements.

In order to gain accurate information on the sample (e.g., the spatial distribution or the diffusion constant of the fluorophore molecules), one has to know as precisely as possible where the fluorescent photons originate from. In other words, one needs to be aware of the detailed spatial profile of fluorescence. This needs the consideration of the following factors:

  • The spatial profile of the excitation light
  • The spatial distribution of the fluorophore molecules also considering diffusion
  • The state transitions (including photobleaching) of the fluorophore molecules according to their photophysical properties without the assumption of a stationary state, but also taking pulse saturation and ground state depletion into account.

In this paper, we present a numerical model which combines all these considerations contrary to previous works that combined some but not all of the mentioned elements. The models described in the preceding publications can be divided into two groups: models in the first group focus on state transitions but without taking into account the spatial profile,911 and some of them also lack the consideration of photobleaching.12 Models in the second group deal with the spatial profile; however, they do not distinguish real photophysical states, but only an “active” and a photobleached state (and in one case, also triplet dynamics), which corresponds to the case in which photobleaching affects molecules in the ground state.2,3,1315

Based on our model, we implemented a simulator program which enables the investigation of the spatial distribution of fluorescence and photobleaching in two-photon excited samples. The model arrangement is as follows: a diffusible, initially homogenous solution of a fluorophore is excited by a focused laser pulse train. The laser light evokes two-photon excitation and other state transitions according to the applied photodynamic model.

A substantial difference between the simulations based on the presented model and real experiments is that the simulator uses a continuum model which determines the probabilities and expected values of the events instead of handling the discrete stochastic events (e.g., photon emission and photobleaching) that occur in practice. Moreover, the simulator can reveal the spatial distribution of photon emission within the excitation volume unlike a real two-photon microscope, which counts the photons emitted from the whole excitation region without considering their precise source location.

The simulated arrangement can serve as a basic model for several measurement and imaging techniques such as the aforementioned two-photon and super-resolution microscopy, FCS, and FRAP.

In Sec. 2, we present our computational model in detail: the model of the illumination, the photodynamic model of the fluorophore, and later on, the diffusion model and the model of the photon detection. We also recall the underlying physical principles. Afterward, in Sec. 3, we present our simulations. Finally, we discuss the results (Sec. 4) and close the paper with conclusions and acknowledgments.

We modeled the following arrangement: a pulsed laser beam is focused on a homogenous fluorophore solution in which excitation, relaxation processes, photon emission, and photobleaching take place (Fig. 1). The sample is divided into cubic cells among which diffusion occur. The simulator calculates the number of fluorophore molecules in each photodynamic state (see Sec. 2.2.2) and the number of emitted fluorescent photons for every volume cell in each time step.

Graphic Jump LocationF1 :

Scheme of the simulation arrangement. The red hour glass shape denotes the illumination laser beam. The dotted rectangular box shows a possible simulation volume instead of which a cylindrical volume is used: the actually simulated volume cells are denoted by the small blue cubes. Axis z points toward the direction of light propagation.

The arrangement has a cylindrical symmetry with respect to the optical axis of the illuminating laser beam. Thus, a volume cell can be referred to by two cylindrical coordinates, namely the radial coordinate r denoting the distance from the beam axis and the axial coordinate z denoting the distance measured from the focus along the beam axis; the origin of the coordinate system is the focus of the illuminating beam. Taking advantage of the cylindrical symmetry, computational complexity can be remarkably reduced: there is no need to simulate the whole illuminated region (denoted by the dotted rectangular box in Fig. 1), but it is sufficient to simulate only one radial slice of it,2 which reduces the number of volume cells from O(n3) to O(n2). The arrangement is also mirror symmetric with respect to the focal plane z=0, which enables a further reduction of the number of volume cells by a factor of 1/2.

Model of the Illumination

A monochromatic, continuous wave laser beam is commonly described as a Gaussian beam (also called Gaussian–Lorentzian profile,9 which is a solution of the paraxial Helmholtz equation).16 However, to generate a photon flux that is intense enough to provoke two-photon excitation, mode-locked lasers are applied, which emit pulses with a duration in the femtosecond range. For instance, titanium:sapphire lasers are prevalently used in two-photon microscopes. They emit 100-fs-long laser pulses with a wavelength between 700 and 1050 nm,6 which corresponds to about 30 to 40 cycles of the electromagnetic wave. Assuming that the temporal shape of these pulses is a Gaussian function, thus their Fourier transform is also a Gaussian function,17 and we can determine that the full-width at half-maximum (FWHM) bandwidth of such a Gaussian pulse is about 55 THz, which corresponds to 15% of the typical central frequency of 375 THz (800 nm). Consequently, the monochromatic continuous wave Gaussian beam offers an acceptable approximate description for such pulses.

In this work, we modeled the illuminating light as a train of focused laser pulses, whose temporal profile is Gaussian (Fig. 2): Display Formula

P(t)=Pmaxexp(t22σpulse2).(1)
In which Display Formula
σpulse=τpulse22ln2(2)
denotes the standard deviation of the Gaussian pulse with FWHM of τpulse. Pmax is the peak power of the pulses, which can be determined from the long-time-averaged laser power Pavg as Display Formula
Pmax=Pavg2ln2πτrepτpulse,(3)
where τrep is the pulse repetition time (i.e., the inverse of the pulse repetition rate).

During the simulations, we consider the laser beam to be monochromatic. (However, as pulses get shorter, the error of this approximation increases, the spectrum of the pulses expands, and the space–time profile of the beam varies significantly. For the description of such pulsed beams, see the work of April.18)

In the Gaussian beam, the light intensity at a point (r,z) is given by the equation Display Formula

I(r,z)=I0[WE0WE(z)]2exp[2r2WE2(z)],(4)
where Display Formula
WE(z)=WE01+(zzE0)2,(5)
in which Display Formula
zE0=πWE02λexc(6)
is the Rayleigh range, Display Formula
WE0=λexcπθ,(7)
is the waist radius, λexc is the wavelength of the excitation light, and θ is the beam divergence.16 Subscript E refers to excitation for the sake of the distinction from the detection profile (see Sec. 2.4). For a continuous wave Gaussian beam, the light intensity at the focus can be determined from the laser power P as Display Formula
I0=2PπWE02.(8)

We suppose that this relationship also remains approximately valid for our pulsed case, i.e., we can assume that in Eqs. (4) and (8), the instantaneous intensity I0(t) at the focus as well as the instantaneous intensity I(r,z,t) at point (r,z) follows the laser power P(t) without delay. Therefore, the light intensity at a point (r,z) and time instance t can be calculated as Display Formula

I(r,z,t)=I0(t)[WE0WE(z)]2exp[2r2WE2(z)],(9)
moreover Display Formula
I0(t)=2P(t)πWE02.(10)

In our model, we neglect the scattering and the attenuation of light in the sample. By the use of the Gaussian beam, the illumination is described as a scalar field of light intensity, i.e., only the magnitude of the electric field is considered, but not its direction. We assume that the orientation of the fluorophore molecules is uniformly random, thus the influence of the orientation is involved in the absorption cross-section. Furthermore, we assume that the lifetime of the relaxation process coupled to photon emission is much longer than the Brownian rotational lifetime of the fluorophore molecules; therefore, the emitted fluorescent light is considered to be depolarized.19

Photodynamic Model of the Fluorophore Molecules
Theoretical description of the rates of one- and two-photon absorptions

In the following, we recall the theoretical background of one- and two-photon absorptions according to the books of Boyd20 and Masters and So1 using the semiclassical model for the description of light–matter interaction. In this approach, the matter is quantized (it is built of particles), whereas light is handled as an electromagnetic field. In the absence of radiation, one can describe the atomic system with the Schrödinger equation: Display Formula

[p22m+V(r)]Ψ(r,t)=iΨ(r,t)t,(11)
whose eigenfunctions and eigenvalues give the possible states and energy levels of the system, where Ψ is the state function, p is the momentum, m is the mass of the particle, V(r) is the potential energy at point r, i is the imaginary unit, and is the reduced Planck constant. According to the electric dipole approximation, illumination of the system can be interpreted as a perturbation during which the electric field induces electric dipoles. Consequently, a new term is to be added to the Schrödinger equation when the atomic system is illuminated: Display Formula
[p22m+V(r)ejE(ω,t)·rj]Ψ(r,t)=iΨ(r,t)t,(12)
where rj is a vector connecting the separated charge ej in the j’th electric dipole induced by electric field E(ω,t) of angular frequency ω at time t. Assuming that the illuminating light is monochromatic and linearly polarized, from Eq. (12) applying the first-order time-dependent perturbation theory, one can express the rate of the transition from the ground state m to an excited state n due to one-photon absorption as Display Formula
Rmn(1)=σmn(1)(ω)I,(13)
where σmn(1)(ω) is the one-photon absorption cross-section in units of cm2/photon, and I is the light intensity or, more precisely, photon flux, measured in units of photon/(cm2s). Similarly, in the case of two-photon absorption, the second-order time-dependent perturbation theory results in the following rate for the transition from the ground state m to an excited state n, via a virtual state k: Display Formula
Rmn(2)=σmn(2)(ω)I2,(14)
where σmn(2)(ω) is the two-photon absorption cross-section in units of cm4s/photon2. The absorption cross-sections can be measured or they can be determined theoretically by Fermi’s golden rule: Display Formula
σmn(1)(ω)=4π2ωnc|μmn|ρ(ω=ωmn)(15)
and Display Formula
σmn(2)(ω)=8π3ω2n2c2|kμmkμknωknω|2ρ(2ω=ωmn),(16)
where ω denotes the angular frequency of the illuminating light, n is the refractive index, c is the speed of light, and ρ is the density of states. The theoretical derivation of the absorption cross-sections requires the knowledge of the transition moment μmn, which can also be calculated theoretically, namely from the wave functions φm and φn of the two states: Display Formula
μmn=jejφm*(0)rjφn(0)d3r.(17)

Nevertheless, for larger molecules, the calculation of the state functions is practically impossible because of the computational complexity of the problem. In our simulations, the cross-sections are to be given as input parameters.

Finally, we note that the connection between the light intensity and the electric field is Display Formula

I=nc2πω|E|2.(18)

Photodynamic states and transitions in our model

We model the states and state transitions by a multilevel photodynamic model taken from the literature911 with modifications. The model includes four different photobleaching routes starting from the excited states. It contains the following states (Fig. 3): ground state S0, excited state S1, triplet state T1, and four photobleached states B1, B2, B3, and B4. We neglect vibrational and rotational states; therefore, we consider the electronic states to be discrete (sharp). In the model, the following transitions are possible:

  • S0S1, S1B2, and T1B4 via two-photon absorption characterized by absorption cross-sections δexc, δpb2, and δpb4,
  • S1B1 and T1B3 via one-photon absorption characterized by absorption cross-sections σpb1 and σpb3,
  • S1S0 radiative relaxation followed by emission of a fluorescent photon, characterized by time constant τf,
  • S1S0 nonradiative relaxation (internal conversion) characterized by time constant τIC,
  • S1T1 intersystem crossing characterized by time constant τISC, and
  • T1S0 nonradiative relaxation characterized by time constant τTS.

Graphic Jump LocationF2 :

The Jablonski diagram of the fluorophore represents the photodynamic states and the possible state transitions of the model. S0, S1, and T1 denote the ground state, excited state, and triplet state, respectively; B1, B2, B3, and B4 indicate the photobleached states. Solid red arrows pointing upward denote photon-absorption-induced transitions (group 1), arrows pointing downward denote the relaxation transitions (group 2) with (solid blue line) or without photon emission (dashed black lines), and k denotes the rate constant of each transition.

We omit stimulated emission as it is neglected also in several referenced papers.9,11,19 The transition from state A to state B is described as follows: the time derivative of the number of molecules in state A in volume cell (r,z) at time instant t (supposing the change is due only to transition AB) is Display Formula

A˙(r,z,t)=kAB(r,z,t)A(r,z,t),(19)
whereas the time derivative of the number of molecules in state B (supposing the change is due only to transition AB) is Display Formula
B˙(r,z,t)=kAB(r,z,t)A(r,z,t),(20)
where kAB is the rate constant of the transition. If more than one transition affect a state, then the changes of the numbers of molecules deriving from each transition are added up.

The photodynamic processes in the model can be divided into two groups:

  1. Transitions in the first group are induced by photon absorption. The rate constant of such a transition depends on the illumination intensity either linearly (kpb1 and kpb3) or quadratically (kexc, kpb2, and kpb4), and it is accordingly proportional to either the one- (σpb1 and σpb3) or two-photon absorption cross-section (δexc, δpb2, and δpb4) which characterizes the process. The rate constant can be calculated as Display Formula
    kAB(r,z,t)=σABλexchcI(r,z,t),(21)
    if the transition is induced by one-photon absorption, and as Display Formula
    kAB(r,z,t)=12δAB(λexchc)2I2(r,z,t),(22)
    if it is coupled to two-photon absorption,11 where σAB and δAB are the one- and two-photon absorption cross-sections at wavelength λexc, respectively, h is Planck’s constant, c is the velocity of light, whereas I(r,z,t) is the illuminating light intensity at point (r,z) at time instance t in units of W/cm2. The space- and time-dependences of the rate constants arise from the space- and time-dependence of the illumination intensity. In Eq. (22), the factor 1/2 expresses that two photons are required for the two-photon excitation process.
  2. The rest of the transitions (relaxation transitions) are intensity-independent and thus space- and time-independent. Each of them is characterized by a global time constant. The rate constant of such a process is determined as the inverse of the time constant of the process (kf=1/τf, kIC=1/τIC, kISC=1/τISC, and kTS=1/τTS).

We suppose that the pulse duration is much shorter than the time constant of the relaxation processes (i.e., τpulseτIC, τf, τISC, τTS), furthermore the pulse length is much shorter than the pulse repetition time (i.e., τpulseτrep). Therefore, during a light pulse, the relaxation processes can be neglected and it can be assumed that a fluorophore molecule is excited atmost once per pulse.911,19 Consequently, it is reasonable to construct two differential equation systems: the first one describes the photon-absorption-induced processes taking place during the light pulses (corresponding to group 1); the second one expresses the relaxation transitions occurring during the interpulse periods (corresponding to group 2).10 Thus, the following two differential equation systems have to be constructed for each volume element: Display Formula

[S˙0S˙1T˙1B˙1B˙2B˙3B˙4]=[kexc00kexc(kpb1+kpb2)000(kpb3+kpb4)0kpb100kpb2000kpb300kpb4]·[S0S1T1](23)
and Display Formula
[S˙0S˙1T˙1]=[(kIC+kf)kTS(kIC+kf+kISC)0kISCkTS]·[S1T1].(24)

The vectors on the right side contain the number of molecules in the given cell in the corresponding states; the vectors on the left are built from the time derivatives of the molecule numbers; whereas the matrices contain the rate constants characterizing the transitions. The molecule numbers and their time derivatives in both equations, and the rate constants in Eq. (23), depend on spatial coordinates (r,z) and time t; however, for the sake of simplicity, it is not signed here [as opposed to Eqs. (19) and (20)]. On the other hand, rate constants in Eq. (24) are independent of space and time as already mentioned.

Computation of the transition probabilities

Based on the ordinary differential equation system formulated in Eq. (23), one can calculate the probability of each photon-absorption-induced state transition during one laser pulse for every volume cell. Because the pulse train consists of identical laser pulses, it is enough to calculate these probabilities once. S0, S1, and T1 are the three states from which transitions are possible. Accordingly, the simulator solves the differential equation system with three different initial conditions formulating that the molecule starts from state S0, S1, or T1. In our model, the time profile of the illumination intensity is Gaussian; thus, it has an infinite support. Nevertheless, the numerical integration is performed on the arbitrarily chosen domain spanning from 4σpulse to +4σpulse, where σpulse denotes the standard deviation of the Gaussian function describing the laser pulse.

The relaxation transitions [together with the diffusion (see Sec. 2.3)] taking place between two consecutive laser pulses are simulated in N steps; that is, the simulation time step is τstep=τrep/N. Similar to the photon-absorption-induced transitions, the simulator calculates the probability of each relaxation transition occurring between laser pulses by solving Eq. (24) numerically for t=τstep. The solution is performed with two different initial conditions which describe that the molecule starts from state S1 or T1. [For Eq. (24), there exists also an analytical solution, namely exponential functions. Transition S1T1S0 is a two-step decay, for which Bateman gave a solution.21]

The numerical solution of the differential equation systems is performed by the ode45 MATLAB function applying the Dormand-Prince method, which combines explicit fourth and fifth order Runge-Kutta formulae.22

As a result, the simulation of a state transition is simplified to the transfer of a certain fraction of the molecules from one state to another according to the probability calculated initially. Pulse saturation and ground state depletion are inherently taken into consideration by the model.

Diffusion Model

As explained above, the sample is modeled as an initially homogenous, diffusible fluorophore solution, which is divided into cubic cells. Hereby, we give the detailed model of the isotropic diffusion taking place in the sample.

Theoretical description of diffusion: Fick’s law

The transportation of molecules due to the concentration gradient is described by Fick’s law of diffusion: Display Formula

J=Dgradρ,(25)
where J is the vector of particle flux density, D is the diffusion constant, and ρ is the concentration of the particles.23 The relation between the diffusion flux and the concentration is given by the continuity equation, which expresses the conservation of the particles: Display Formula
divJ=ρt.(26)

From Eqs. (25) and (26), the diffusion equation can be obtained: Display Formula

ρt=DΔρ.(27)

The solution of Eq. (27) is Display Formula

ρ(x,y,z,t)=1(4πDt)3exp(x2+y2+z24Dt),(28)
which describes the concentration evolved by diffusion at point (x,y,z) at time t if the matter is placed at the origin, t=0.23 Separating the diffusion along the three axes, we obtain Display Formula
ρ(x,y,z,t)=14πDtexp(x24Dt)·14πDtexp(y24Dt)·14πDtexp(z24Dt),(29)
i.e., the concentration of the molecules along each axis is described by a Gaussian distribution with zero mean and variance of 2Dt.

Discretized diffusion model

Simulation needs discretization both in space and time. The spatial discretization is achieved by the cubic mesh described at the beginning of Sec. 2, whereas the time step of the diffusion simulation is τstep=τrep/N as was stated in Sec. 2.2.3. (In Sec. 3.1.3, we present that even N=1 can be a convenient choice.) In every time step, a molecule can stay at the same volume cell or can diffuse to one of the two neighboring cells along each axis; the probability of diffusion is the same toward the negative and positive directions. Let X be a discrete random variable which denotes the dislocation of a molecule from the origin along axis x. It can take three values:

  • x1=Δx: the dislocation of the molecule is one volume cell toward the negative direction along the axis,
  • x2=0: the molecule stays at the same volume cell, and
  • x3=Δx: the dislocation of the molecule is one volume cell toward the positive direction along the axis.

The probabilities that belong to these three events are the following:

  • p1=Pr(X=x1)=Pr(X=Δx)=p,
  • p2=Pr(X=x2)=Pr(X=0)=12p, and
  • p3=Pr(X=x3)=Pr(X=Δx)=p,

where 0<p<(1/3). Therefore, the expectation value and variance of X are Display Formula

E[X]=i=13pixi=0(30)
and Display Formula
Var[X]=E[(XE[X])2]=i=13pixi2=2p·(Δx)2.(31)

Let Y(k) be a random variable which denotes the dislocation of a molecule after k time steps, i.e., Display Formula

Y(k)=X+X++Xk=k·X.(32)

According to the central limit theorem,17 as k, Y(k) is asymptotically normal (Gaussian) with mean Display Formula

E[Y(k)]=μ=i=1kE[X]=k·E[X]=0(33)
and variance Display Formula
Var[Y(k)]=σ2=i=1kVar[X]=k·Var[X]=2k·p·(Δx)2.(34)

Equation (29) states that the distribution of the diffusing matter is Gaussian with zero mean and a variance of 2Dt; thus, for the variance σ2 of Y(k), equation Display Formula

σ2=2k·p·(Δx)2=2Dt,(35)
should hold. Substituting t=kτstep into Eq. (35), we obtain Display Formula
p=DτstepΔx2,(36)
for the probability that a molecule moves to one of the neighboring volume cells along a coordinate axis during one time step.

Cylindrical symmetry

The three axes along which we describe the diffusion are axes r, z, and a third axis y, which is perpendicular to both axes r and z. Therefore, in order to simulate the diffusion along axis y, we need two adjacent volume cell layers lying next to the simulation layer (Fig. 4). Because of the cylindrical symmetry with respect to axis z, these adjacent layers are identical. The number of molecules in the cells of these layers is estimated by interpolation based on the values of the simulation layer. For the sake of computational simplicity, linear interpolation is used.

Graphic Jump LocationF3 :

The estimation of the neighboring cell layer for the calculation of diffusion is achieved by linear interpolation using the values of the simulation layer. Axis z points from the reader into the figure at point O.

Because of the cylindrical symmetry, the concentration at point P has to be equal to the concentration at point Q since both points are at distance r0 from the origin. From the concentrations cA and cB at points A and B, the concentration cQ=cP at points Q and P is determined by linear interpolation: Display Formula

cP=cQcA+(r0r0)(cBcA),(37)
where denotes the floor function. Because the distance between points A and B is unity, r0 is equal to the distance between points O and A.

Boundary conditions

During the simulation of diffusion, we use zero-flux (Neumann) boundary condition at the outer edges, i.e., the edge layer is repeated and is put beside the simulation region. The cells situated on axes r and z form “virtual” edges because they are actually in the middle of the simulated region. Therefore, in these cases to preserve the delineated symmetry of the arrangement, not the edge layer (which is located on the axis of symmetry), but the next one is repeated and put beside the simulation region during the simulation.

Model of the Photon Detection

Here, we present a possible model for photon detection, although in this paper we do not cover the detected number of photons nor the detection profile.

The probability that one detects a fluorescent photon originating from the volume element located at (r,z) is Display Formula

pdet(r,z)=η·Φ·Y(r,z),(38)
where η is the efficiency of the detection, Φ is the fractional solid angle of the observation, and Y(r,z) is the observation beam profile.19 In the case of full aperture detection, Display Formula
Φ=NA24n2(39)
and Display Formula
Y(r,z)=1.(40)

It means that the fluorescent photons are detected with the same probability, independently from the point from where they were emitted. (NA denotes the numerical aperture, and n is the refractive index of the sample). However, for confocal observation through a single-mode fiber, Display Formula

Φ=λemit24π2n2WD02(41)
and Display Formula
Y(r,z)=[WD0WD(z)]2exp[2r2WD2(z)];(42)
i.e., the observation beam profile is described by a Gaussian beam, whose waist radius is WD0. λemit denotes the wavelength of the emitted fluorescent light which is to be detected. (Subscript D refers to detection.).

Process of the Simulation

In the beginning of the simulation, the probability of each state transition is calculated for every volume cell (see Sec. 2.2.3). Then comes the simulation of the state transitions and photon emission evoked by the pulsed illumination. Initially, all the molecules are in the ground state (S0). At the end of the first light pulse, a fraction of the molecules given by the transition probabilities is transferred to other states by the photon-absorption-induced transitions (group 1 in Sec. 2.2.2). Until the start of the second light pulse, relaxation processes (group 2 in Sect. 2.2.2) and diffusion take place. They are simulated in the same way: some fractions of the molecules determined by the transition probabilities are transferred to other states or to the neighboring cells. These two phases compose one laser cycle and they are repeated for the required times during the simulation process (Fig. 2). The pseudocode of the algorithm is provided in Appendix.

Graphic Jump LocationF4 :

The temporal shape of the laser pulse train and the main steps of the simulation. Each laser pulse is described as a Gaussian function with a full-width at half-maximum length of τpulse. The time period between two consecutive pulses is τrep. Note that in the figure the ratio of these two periods is distorted for the sake of demonstrativeness; in fact, τpulse is much shorter than τrep. Transitions related to photon absorption taking place during the laser pulses are simulated in one step, whereas the relaxation transitions and diffusion occurring between the laser pulses are simulated in N (1 or more) steps.

Selection of the Simulation Parameters

Our goal here is, illustrating the capabilities of our model and simulator, to demonstrate by characteristic simulation cases how different factors (laser power, photobleaching cross-section, diffusion coefficient, and intersystem crossing) affect fluorescence. A notable advantage of the simulator is that it makes possible the separate investigation of these factors.

We tried to choose the simulation parameters in such a way that they could reflect realistic circumstances. Nevertheless, they do not refer to a given measurement on a real fluorophore by a real device. First, we fixed the following parameters: pulse length τpulse=100fs, pulse repetition time τrep=12.5ns (i.e., the pulse repetition rate was 80  MHz), and the wavelength of the illuminating light λexc=800nm, whose values correspond to the Ti:sapphire laser commonly used in two-photon microscopes. Unless indicated, the diffusion and the relaxation processes are simulated in one step per laser cycle; i.e., N=1 and τstep=τrep. We set θ=60deg for the beam divergence. In every simulation, we used the arbitrarily chosen 1μmol/L for the concentration of the homogenous fluorophore solution; however, this value behaves as a simple multiplier when molecule and photon numbers are calculated by the simulator. We set the fluorophore parameters to be in the same order of magnitude as values found in the literature10,11 for real fluorophores. We used 2×1049cm4s=20GM for the excitation two-photon absorption cross-section δexc and 3 ns for the lifetime τf of fluorescent relaxation. In each presented simulation, we disabled either all the photobleaching routes (by setting their absorption cross-sections to zero) or we enabled one of them (by setting its absorption cross-section to some nonzero value given later). We also disabled internal conversion in all the simulations and intersystem crossing in all but in one virtual experiment setting their lifetimes to τIC= and τISC=. (The one exception is indicated later, see Sec. 3.2.5.)

Calibration of the laser power

First, we ran simulations in order to examine the spatial profile of the excitation probability at different values of the laser power Pavg. In this experiment, we disabled photobleaching (i.e., we set σpb1, δpb2, σpb3, and δpb4 to 0) and diffusion (i.e., we set diffusion constant D to 0). We set the radius and height of the simulated cylinder to R=5μm and H=10μm, respectively, and the cell size to Δx=0.01μm. We used T=1.25×106s for the simulation length, which corresponds to 100 laser cycles.

The probability of excitation at the focus during one laser pulse increased nearly linearly with the laser power in the case of moderate illumination (5–30  mW) [Fig. 5(a)]. At higher laser powers, the excitation probability saturated, but the number of emitted photons (i.e., the signal) still increased [Fig. 5(b)]. Nevertheless, when we increased the illumination intensity, the fluorescing region expanded [Figs. 5(c) and 5(d)]. It means that the uncertainty of the source location of the fluorescent photons undesirably increased. Based on these results, we chose two laser power levels for the further simulations: 25 and 40 mW. They correspond to the cases when the probabilities that a fluorophore molecule at the focus is excited during one pulse are 70.4% and 95.6%, respectively. Figure 6 depicts the spatial profile of the probability that a fluorophore molecule is excited during one laser pulse in the cases of the mentioned two laser power values.

Graphic Jump LocationF5 :

Dependence of (a) the excitation probability at the focus, (b) the number of emitted photons, and (c) and (d) the excitation profile on the illumination laser power. (c) and (d) The probability along axes r and z that a fluorophore molecule is excited during one laser pulse. The following values were used for the laser power: 5, 10, 15, 20, 25, 30, 35, 40, 50, and 60 mW. The higher the laser power is, the darker and the higher the curve is. Note the difference between the scales of the horizontal axes of (c) and (d).

Graphic Jump LocationF6 :

The spatial profile of the probability that a fluorophore molecule is excited during one laser pulse in the cases when the laser power Pavg is (a) 25 mW and (b) 40 mW.

Calibration of the cell size

In the next simulation, we checked to what degree the cell size influences the results. Because the illumination intensity changes faster in the space near the focus, whereas its gradient is smaller on the periphery, we ran simulations of length T=1.25×106s on a small region (R=0.32μm and H=0.96μm) around the focus with different cell sizes. Photobleaching and diffusion were still disabled. Figure 7(a) summarizes the results, whereas Fig. 7(b) depicts the number of simulated cells, which is proportional to the computation time and the memory requirement. Since it was a reasonable trade-off between accuracy and computational time, we retained Δx=0.01μm for the cell size in further simulations. [For this cell size, the error of the emitted-photon number is 0.60% and 0.69% for the two chosen laser power values (25 and 40 mW) considering the case with Δx=1.25×103μm as the reference.]

Graphic Jump LocationF7 :

Dependence of (a) the fluorescent photon number and (b) the number of simulated cells on the cell size. In (a), the reference is the photon number obtained for the simulation case with Δx=1.25×103μm.

Calibration of the simulation of diffusion

Here, we investigate how much it modifies the emission profile if the time step τstep of the simulation of diffusion (and relaxation transitions) is reduced. The influence of the diffusion on the fluorescent profile increases with the diffusion constant, the degree of photobleaching, and the length of the simulation. Accordingly, to study the “worst case,” we used Pavg=25mW and σpb1=5×1022cm2 for the laser power and the photobleaching cross-section, respectively, since, among those presented in this paper, this is the setting that results in the highest degree of photobleaching. Furthermore, we set the diffusion constant to the highest value occurring in the presented simulations, namely D=8×1010m2/s. The simulation length was T=1.25×103s (100,000 cycles), and we set the radius and the height of the simulation volume to R=0.5μm and H=1.5μm, respectively. For N, which denotes in how many steps per laser cycle the diffusion and relaxation transitions are simulated, we used eight different values: 1,2,4,,128. In other words, the simulation time step τstep was set to the pulse repetition time τrep=12.5ns of the mode-locked laser and its 1/2,1/4,,1/128.

Figures 8(a)8(c) delineate the local relative error of the number of emitted photons per volume cell in the 100,000th cycle in the cases when diffusion and relaxation processes were simulated in 1, 2, and 4 steps per laser cycle. Figure 8(d) shows the maximum local relative error of the number of photons emitted in the 100,000th cycle as a function N. In each case, the reference is the photon emission profile obtained for the 128-steps-per-cycle case. These diagrams demonstrate the convergence of the calculated photon emission distribution as the simulation time step τstep tends to zero.

Graphic Jump LocationF8 :

(a)–(c) The local relative error of the number of emitted photons per volume cell in the 100,000th cycle in the cases when diffusion and relaxation processes were simulated in 1, 2, and 4 steps per laser cycle. (d) The maximum local relative error of the number of photons emitted in the 100,000th cycle as a function of the number of simulation steps per laser cycle. In each case, the radius and height of the simulation volume were R=0.5μm and H=1.5μm, respectively, and the reference is the photon emission profile obtained for the 128-steps-per-cycle case. Note the difference between the scales of (a)–(c).

Even in the one-step-per-cycle case [Fig. 8(a)], the relative error was less than 1‰ in the majority of the simulation region. There were two exceptions: some of the cells located at the boundary of the simulation volume (here, the relative error was about 2%) and the outer regions located near the focal plane (where the relative error was about approximately 1%). However, the fact that the rate of excitation (and thus that of photon emission and photobleaching) is small in these regions (see Fig. 6) diminishes the significance of this error. Consequently, even in the presented “worst case” (high photobleaching and intense diffusion), it is enough to simulate the diffusion in one step per laser cycle. Therefore, we used this setting in the following simulations.

Now, let us consider the error deriving from the finiteness of the simulation volume. The applied zero-flux boundary condition works properly if the concentration gradients are small at the edge regions. Since the spatial inhomogeneity of the molecules in different states originates from the photophysical transitions taking place mostly at the focus, it can be assumed that increasing the size of the simulation volume decreases the error deriving from the finiteness of the computation region. Thus, to find the appropriate dimensions, we ran simulations in which the radius R of the simulation volume was swept from 0.5μm to 3.0μm, while the height H of the simulation volume was, in each case, three times larger than the radius.

Figures 9(a)9(c) depict the local relative error of the number of emitted photons per volume cell in the 100,000th cycle in the cases of the three smallest simulation volumes. Figure 9(d) shows the maximum local relative error of the number of photons emitted in the 100,000th cycle within the inner region (of radius R=0.5μm and height H=1.5μm) as a function of the radius of the simulation volume. In each case, the reference is the photon emission profile obtained for the largest simulated volume (of radius R=3μm and height H=9μm). These diagrams confirm that as the simulation volume is expanded, the maximum of the local relative error falls rapidly. When the radius and height of the simulation volume are set to R=1.5μm and 4.5μm, the error is tolerably small (approximately 2‰ to 3‰); thus, this setting can be regarded as an acceptable trade-off between accuracy and simulation time. If the diffusion constant is smaller, then the degree of the error diminishes to a great extent. It would, therefore, be reasonable to choose a smaller simulation volume in these cases. Nonetheless, for the sake of simplicity and consistency, we set the dimensions of the simulation volume to the aforementioned values (R=1.5μm and 4.5μm) in all those simulation cases when diffusion is enabled.

Graphic Jump LocationF9 :

(a)–(c) The local relative error of the number of emitted photons per volume cell in the 100,000th cycle in the cases when the radius R of the simulation volume was 0.5μm, 1.0μm, and 1.5μm, respectively, and the height H of the simulation volume was three times its radius. (d) The maximum local relative error of the number of photons emitted in the 100,000th cycle within the inner region (of radius R=0.5μm and height H=1.5μm) as a function of the radius of the simulation volume. In each case, the reference is the photon emission profile obtained for the largest simulated volume (of radius R=3μm and height H=9μm). Note the difference between the scales of (a)–(c).

Investigation of the Three-Dimensional Fluorescent Profile and Its Time Evolution
Factors affecting the photon emission profile

We have arrived at the simulations investigating the spatial profile of fluorescence and its changes over time due to photobleaching. Figure 10 depicts 12 simulation cases in which we set different input parameters. We used two different values for laser power Pavg (25 and 40 mW), diffusion constant D (zero and 2×1010m2/s), and for photobleaching cross-section σpb1: 5×1023cm2 (a value found in the paper of Niesner et al.10) and a value 10 times larger (5×1022cm2) to simulate increased photobleaching. (The three other photobleaching routes were disabled, i.e., their absorption cross-sections were set to 0.) The heat maps show the normalized number of emitted photons per volume cell during the 10th, 10,000th, and 100,000th laser cycles: in each row, the maximal local value of the fluorescent photon number corresponds to 1.

Graphic Jump LocationF10 :

Spatial distribution of photon emission. The images depict the normalized number of emitted photons per volume cell during the 10th, 10,000th, and 100,000th laser cycles. The normalization was carried out row by row: in each row, the maximal local value of the photon number corresponds to 1. In the middle row [(e)–(h)], both the laser power and the photobleaching cross-section σpb1 were moderate. Even in this case, the decline of fluorescence due to photobleaching was perceptible after an illumination of 10,000 laser cycles (f), and it became severe after 100,000 cycles (g), especially at the focus, where a “dark hole” evolved. Laser power was increased in the first row [(a)–(d)] and the photobleaching cross-section in the third one [(i)–(l)]. In both cases, photobleaching was accelerated. In columns 1–3, diffusion was disabled; however, it occurred in column 4 compensating photobleaching more or less. Note that in the simulations that are presented in column 4 and in the case for which diffusion was enabled, the radius and height of the simulated region were R=1.5μm and H=4.5μm, respectively; however, only the focal part of the profile is shown here.

Even when both the laser power and photobleaching cross-section σpb1 were moderate (25 mW and 5×1023cm2, middle row of the table), the decline of fluorescence due to photobleaching was perceptible after a 1.25×104-s-long illumination [10,000 laser cycles, Fig. 10(f)], and it became severe after 1.25×103s [100,000 cycles, Fig. 10(g)], especially at the focus, where a “dark hole” evolved. When the photobleaching cross-section was set to 10 times higher (third row of the table), then after 10,000 cycles, we got a profile [Fig. 10(j)] which is similar to the one obtained in the case of moderate photobleaching cross-section after 100,000 cycles [Fig. 10(g)]. Not surprisingly, the increase of the laser power also accelerated photobleaching (first row of the table). Diffusion, however, compensated photobleaching more or less (compare columns 3 and 4). (This phenomenon is the foundation of FRAP.)

Figure 11 shows the emission profiles depicted in Figs. 10(e) and 10(g) in a three-dimensional plot.

Graphic Jump LocationF11 :

Emission profile in the (a) 10th and (b) 100,000th laser cycles. The results depicted here are the same as those in Figs. 10(e) and 10(g). Note the difference between the scales of the vertical axes of (a) and (b).

Time evolution of the fluorescent profile in the cases of one- and two-photon-induced photobleaching

In the following simulations, we investigated the time evolution of the spatial profile of fluorescence in the cases of one- and two-photon-induced photobleaching occurring via state S1. For the former transition, we used σpb1=5×1023cm2 again, whereas we set δpb2 to 6.534×1054cm4s in order that the probability that a fluorophore molecule being in state S0 and located at the focus is photobleached during one laser pulse would be approximately the same for the two cases. In these simulations, we set the laser power to Pavg=40mW and disabled the diffusion.

The results are summarized in Fig. 12. The first row of the table [Figs. 12(a)12(d)] depicts the number of emitted photons along axes r and z during the 1st, 2500th, 5000th, 7500th, 10,000th, 20,000th, 40,000th, and 100,000th laser cycles. The darker the color of the curve, the larger the number of the preceding laser cycles. As more and more laser pulses hit the sample, the number of fluorescent photons falls: we can observe the evolution of the “dark hole” at the focus [see also Figs. 10(c), 10(g), 10(j), and 10(k)]. Another noticeable feature is the appearance of the small “hills” in the fluorescence profile next to the focus. Their width is much smaller than that of the initial fluorescent profile depicted by the uppermost grayish curve in each plot. A distinct difference between the one- and two-photon cases is that two-photon-induced photobleaching is restricted to a smaller volume; therefore, the “walls of the hole” are steeper in this case.

Graphic Jump LocationF12 :

Spatial distribution of fluorescence and its time evolution caused by one- and two-photon-induced photobleaching. (a)–(d) The number of emitted photons along axes r and z during the 1st, 2500th, 5000th, 7500th, 10,000th, 20,000th, 40,000th, and 100,000th laser cycles. The darker the color of the curve, the latter the state is. (e) and (f) The number of emitted photons as a function of the distance of the photon source from the focus. (g) and (h) The characteristic distances of fluorescence over time: r50, r75, and r90 denote the minimal radii of spheres which are located at the focus and contain the volume cells that emit 50%, 75%, and 90% of the fluorescent photons within the sphere of radius R=1.5μm, located at the focus. The average distance of the location of photon emission from the focus is denoted by linked black squares. Note the difference between the scales of the horizontal axes of (a) and (b), as well as (c) and (d).

Figures 12(e) and 12(f) show the number of emitted photons as a function of the distance of the photon source from the focus. Figures 12(g) and 12(h) show the characteristic distances in the case of different pulse numbers (i.e., different illumination durations): r50, r75, and r90 denote the minimal radii of spheres, which are located at the focus and contain the volume cells that emit 50%, 75%, and 90% of the fluorescent photons within the sphere of radius R=1.5μm, located at the focus. Along with the average distance of the location of the photon emission from the focus (ravg), these values can be used to characterize the spatial resolution of the fluorescent method. The trends mentioned in the previous paragraph can also be observed in these diagrams: the number of fluorescent photons (i.e., the signal) decreased gradually in time, especially in the focal region. The average distance of the photon emission from the focus shifted toward the periphery and the characteristic distances increased; therefore, the resolution deteriorated because of photobleaching. This effect was more expressed in the case of the one-photon-induced photobleaching.

Effect of the photobleaching cross-section

Next, we investigated the dependence of the fluorescent photon number on the photobleaching absorption cross-sections. We separately examined the cases of one- and two-photon-induced photobleaching. We used σpb1=5×1023cm2 and δpb2=6.534×1054cm4s as references, and we swept these parameters in the range of three orders of magnitude. The results are summarized in Fig. 13.

Graphic Jump LocationF13 :

Dependence of the number of emitted photons on the absorption cross-sections that describe photobleaching. The photon numbers are depicted relatively to the reference photon numbers obtained by simulations in which the photobleaching cross-section was set to (a) σpb1=5×1023cm2 and (b) δpb2=6.534×1054cm4s (denoted by larger squares). (The radius and height of the simulated region were R=1.5μm and H=4.5μm, respectively.)

When we decreased σpb1 by a factor of 32, the number of emitted photons increased by only 28%, whereas when we increased σpb1 by the same factor, the number of emitted photons fell to 50% of the reference.

We obtained similar results for the case of two-photon-induced photobleaching. In this latter case, the emitted-photon number changed less (namely in the range between 38% and +17% with respect to the reference) compared with the one-photon case, which can be explained by the fact that the two-photon-induced photobleaching is more confined in space than the one-photon-induced one (see Sec. 3.2.2).

Figure 14 shows the decay of photon emission in time due to photobleaching. The thick curves belong to the aforementioned reference cases (σpb1=5×1023cm2 and δpb2=6.534×1054cm4s); the reference photobleaching cross-section values were both doubled (darker curves) and halved (lighter curves) five times.

Graphic Jump LocationF14 :

Time evolution of the emitted-photon number in the case of different photobleaching cross-sections. The thick curve in each diagram belongs to the simulation case in which the photobleaching cross-section was set to (a) σpb1=5×1023cm2 and (b) δpb2=6.534×1054cm4s. These values were both doubled and halved five times: the higher the photobleaching cross-section is, the lower and darker the curve is. (The radius and height of the simulated region were R=1.5μm and H=4.5μm, respectively.)

Effect of the diffusion

To investigate the effect of diffusion on the number of fluorescent photons, we swept the diffusion constant D from 0 to 8×1010m2/s. (For the sake of comparison, the diffusion constant of glycine in water is about 109m2/s.24) Figure 15 shows the dependence of the number of emitted photons on the diffusion constant. The curves with full and empty triangle marks show the photon numbers obtained for the cases when photobleaching did not occur in the sample; that is, these constant curves indicate the maximal photon number for the given illuminations. We used two different laser powers: 25 and 40 mW. In the 25-mW case, if the photobleaching cross-section was moderate (σpb1=5×1023cm2), diffusion was able to compensate photobleaching. But when photobleaching was more intense (in the case of higher laser power and/or higher photobleaching absorption cross-section), then the number of fluorescent photons decreased even if the diffusion coefficient was large. Additionally, Fig. 15 illustrates that the number of fluorescent photons is more sensitive to the change of the diffusion constant if photobleaching is intensified.

Graphic Jump LocationF15 :

Dependence of the number of emitted fluorescent photons on the diffusion constant. (The radius and height of the simulated region were R=1.5μm and H=4.5μm, respectively.)

Effect of the intersystem crossing

In the simulations presented up to this point, we disabled intersystem crossing by setting its lifetime to . Here, we discuss the effect of the opening of this relaxation route. Intersystem crossing means a roundabout way for fluorophore molecules to return from the excited state to the ground state, namely without photon emission; thus, it decreases the quantum efficiency.

The simulation results (Fig. 16) corresponded to the expectations. The smaller the lifetime of intersystem crossing (τISC) is, the fewer photons were emitted, especially if transition T1S0 was slow, i.e., τTS was large. The return from the triplet state T1 to the ground state S0 takes time, hence with every laser pulse, the number of fluorophore molecules in the triplet state increased toward an equilibrium value [Fig. 17(c)]. Meanwhile, also approaching an equilibrium, the number of emitted photons fell, particularly at the focus [Figs. 17(a) and 17(b)]. As a result, the spatial profile of photon emission became flattened.

Graphic Jump LocationF16 :

Dependence of the number of emitted photons (a) on the lifetime of intersystem crossing (τISC) and (b) on the lifetime of the transition T1S0 (τTS). (The radius and height of the simulated region were R=1.5μm and H=4.5μm, respectively.)

Graphic Jump LocationF17 :

Time evolution of fluorescence in the occurrence of intersystem crossing. (a) and (b) The number of emitted photons along axes r and z during the 1st,50th,100th,,2000th laser cycles. The later moment, the darker the color of the curve. (c) The number of emitted photons per laser cycle and the number of molecules in state T1 over time. (The radius and height of the simulated region were R=1.5μm and H=4.5μm, respectively.)