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

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 in two-photon-excited diffusible samples exposed to photobleaching: a simulation-based study 1

Introduction
Fluorescence is a two-phase process during which photon absorption is followed by photon emission. Its high sensitivity and specificity 1 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 (subfemtoliter 6 ) 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, [9][10][11] 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,[13][14][15] 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.

Simulation Model
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.
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ðn 3 Þ to Oðn 2 Þ. 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): (1) In which denotes the standard deviation of the Gaussian pulse with FWHM of τ pulse . P max is the peak power of the pulses, which can be determined from the long-time-averaged laser power P avg as 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   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.
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 I 0 ð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 moreover 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 2.2 Photodynamic Model of the Fluorophore Molecules In the following, we recall the theoretical background of oneand two-photon absorptions according to the books of Boyd 20 and Masters and So 1 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: 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 pointr, 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: p 2 2 m þ VðrÞ − X e jẼ ðω;tÞ ·r j Ψðr;tÞ ¼ iℏ ∂Ψðr;tÞ ∂t ; (12) wherer j is a vector connecting the separated charge e j in the j'th electric dipole induced by electric fieldẼðω; 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 where σ ð1Þ mn ðωÞ is the one-photon absorption cross-section in units of cm 2 ∕photon, and I is the light intensity or, more precisely, photon flux, measured in units of photon∕ðcm 2 sÞ. Similarly, in the case of two-photon absorption, the secondorder 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: where σ ð2Þ mn ðωÞ is the two-photon absorption cross-section in units of cm 4 s∕photon 2 . The absorption cross-sections can be measured or they can be determined theoretically by Fermi's golden rule: and σ ð2Þ mn ðωÞ ¼ 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: 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

Photodynamic states and transitions in our model
We model the states and state transitions by a multilevel photodynamic model taken from the literature 9-11 with modifications. The model includes four different photobleaching routes starting from the excited states. It contains the following states ( Fig. 3): ground state S 0 , excited state S 1 , triplet state T 1 , and four photobleached states B 1 , B 2 , B 3 , and B 4 . We neglect vibrational and rotational states; therefore, we consider the electronic states to be discrete (sharp). In the model, the following transitions are possible: and T 1 → B 4 via two-photon absorption characterized by absorption cross-sections δ exc , δ pb2 , and δ pb4 , • S 1 → B 1 and T 1 → B 3 via one-photon absorption characterized by absorption cross-sections σ pb1 and σ pb3 , • S 1 → S 0 radiative relaxation followed by emission of a fluorescent photon, characterized by time constant τ f , • S 1 → T 1 intersystem crossing characterized by time constant τ ISC , and 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 A → B) is _ Aðr; z; tÞ ¼ −k AB ð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 A → B) is _ Bðr; z; tÞ ¼ k AB ðr; z; tÞAðr; z; tÞ; (20) where k AB 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 (k pb1 and k pb3 ) or quadratically (k exc , k pb2 , and k pb4 ), 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 if the transition is induced by one-photon absorption, and as if it is coupled to two-photon absorption, 11 where σ AB and δ AB are the one-and two-photon absorption crosssections 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∕cm 2 . The space-and timedependences 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 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. [9][10][11]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: 2 6 6 6 6 6 6 6 6 6 6 6 6 4 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 photonabsorption-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. S 0 , S 1 , and T 1 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 S 0 , S 1 , or T 1 . 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 S 1 or T 1 . [For Eq. (24), there exists also an analytical solution, namely exponential functions. Transition S 1 → T 1 → S 0 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: 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: From Eqs. (25) and (26), the diffusion equation can be obtained: The solution of Eq. (27) is 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 ρðx; y; z; tÞ 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. The probabilities that belong to these three events are the following: where 0 < p < ð1∕3Þ. Therefore, the expectation value and variance of X are and Let Y ðkÞ be a random variable which denotes the dislocation of a molecule after k time steps, i.e., According to the central limit theorem, 17 as k → ∞, Y ðkÞ is asymptotically normal (Gaussian) with mean and variance 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 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. 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 r 0 from the origin. From the concentrations c A and c B at points A and B, the concentration c Q ¼ c P at points Q and P is determined by linear interpolation: where b: : : c denotes the floor function. Because the distance between points A and B is unity, br 0 c 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 p det ð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, and Yðr; zÞ ¼ 1: 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, i.e., the observation beam profile is described by a Gaussian beam, whose waist radius is W D0 . λ 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 (S 0 ). 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.

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 ¼ 100 fs, pulse repetition time τ rep ¼ 12.5 ns (i.e., the pulse repetition rate was 80 MHz), and the wavelength of the illuminating light λ exc ¼ 800 nm, 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 θ ¼ 60 deg 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 literature 10,11 for real fluorophores. We used 2 × 10 −49 cm 4 s ¼ 20 GM 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 P avg . 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 × 10 −6 s 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.

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 × 10 −6 s 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

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 P avg ¼ 25 mW and σ pb1 ¼ 5 × 10 −22 cm 2 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 × 10 −10 m 2 ∕s. The simulation length was T ¼ 1.25 × 10 −3 s (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   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 128steps-per-cycle case. These diagrams demonstrate the convergence of the calculated photon emission distribution as the simulation time step τ step tends to zero.
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 Fig. 8 (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). 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.

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 P avg (25 and 40 mW), diffusion constant D (zero and 2 × 10 −10 m 2 ∕s), and for photobleaching cross-section σ pb1 : 5 × 10 −23 cm 2 (a value found in the paper of Niesner et al. 10 ) and a value 10 times larger (5 × 10 −22 cm 2 ) to simulate increased photobleaching. (The three other photobleaching routes were disabled, i.e., their absorption crosssections 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.
Even when both the laser power and photobleaching crosssection σ pb1 were moderate (25 mW and 5 × 10 −23 cm 2 , middle row of the table), the decline of fluorescence due to photobleaching was perceptible after a 1.25 × 10 −4 -s-long illumination [10,000 laser cycles, Fig. 10(f)], and it became severe after 1.25 × 10 −3 s [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.

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 S 1 . For the former transition, we used σ pb1 ¼ 5 × 10 −23 cm 2 Fig. 9 (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).
again, whereas we set δ pb2 to 6.534 × 10 −54 cm 4 s in order that the probability that a fluorophore molecule being in state S 0 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 P avg ¼ 40 mW and disabled the diffusion. The results are summarized in Fig. 12. The first row of the  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. 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. 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): r 50 , r 75 , and r 90 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 (r avg ), 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 × 10 −23 cm 2 and δ pb2 ¼ 6.534 × 10 −54 cm 4 s as references, and we swept these parameters in the range of three orders of magnitude. The results are summarized in Fig. 13.
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-photoninduced 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 onephoton case, which can be explained by the fact that the twophoton-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 × 10 −23 cm 2 and δ pb2 ¼ 6.534× 10 −54 cm 4 s); the reference photobleaching cross-section values were both doubled (darker curves) and halved (lighter curves) five times.

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 × 10 −10 m 2 ∕s. (For the sake of comparison, the diffusion constant of glycine in water is about 10 −9 m 2 ∕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 × 10 −23 cm 2 ), 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.

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. Fig. 12 Spatial distribution of fluorescence and its time evolution caused by one-and two-photoninduced 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: r 50 , r 75 , and r 90 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). experiments is that our model computes probabilities and expectation values of photon numbers while in reality there are discrete, stochastic events: only whole molecules exist and only whole photons are emitted. Thus, the experimentally measured signals contain fluctuations and noise both in time and space, contrary to the simulation results which are "smooth." Theoretically, if the experiments are repeated many times, then the average of their results should approach to the simulation results.
We implemented the simulator program in MATLAB and ran the simulations on a computer with a 2.4-GHz dual-core processor and 4-GB RAM. The run-time of a simulation case depends on the number of volume cells (which is the function of the cell size and the dimensions of the simulated volume), the number of simulated laser pulses, the number of simulation steps per laser cycle, the number of the enabled state transitions, and whether diffusion is enabled. In terms of these factors, our simulation cases varied from each other to a great extent, which resulted in large differences in their run-times. For instance, the emission profiles in Fig. 10 were calculated within less than a minute if the diffusion was not enabled and the dimensions of the simulation volume were R ¼ 0.5 μm and H ¼ 1.5 μm, but it took 16 min when diffusion was enabled and both the radius and the height of the simulation volume were three times larger. Not surprisingly, increasing the number of simulation steps per laser cycle lengthened the simulation extremely: the 128-steps-percycle case resulted in a simulation time of 2.5 h. Presumably, parallelization on a multicore architecture could significantly accelerate the computation, especially when diffusion is disabled in the simulation.
There are several possible application fields of the model: taking into account the effects described above and using more accurate fluorescent profiles might facilitate the interpretation of the experimental data of several fluorescent techniques. The model might help to determine the diffusion constants in FRAP experiments and to solve the inverse problems of twophoton and super-resolution microscopy (i.e., the reconstruction of the fluorophore distribution from the measurement data) more effectively, achieving a higher image quality. In addition, the model might promote the optimization of the operating parameters (e.g., illumination laser power) of the measurement device in a particular experiment. The photophysical model and its parameters are, however, unknown for many fluorophores, which can limit the applicability of the model for such purposes. Possibly, the simulator can also assist in the determination of these parameters.
In this work, we simulated illumination with a fixed focus; nevertheless, in many fluorescent techniques, the sample is scanned with the laser beam or special illumination patterns are used. Adapting the computational model to such experimental arrangements in order to simulate particular techniques needs further work. In these cases, the cylindrical symmetry vanishes; i.e., the computational time and memory requirements increase to a great extent.
The presented model could be easily adapted for the investigation of fluorescence in one-photon-excited samples. However, in this case, as one-photon excitation is less confined in space than two-photon excitation, the size of the simulated volume should be expanded, which would increase the simulation time and the memory requirement.
We have already reported simulation results in this field in a conference paper. 25 At that stage of our work, we used an equation given by Wang et al. 26 to describe the instantaneous illumination light intensity at a given point ðr; zÞ in a pulsed Gaussian beam. The main drawback of Wang's equation is that for large values of the radial coordinate r, the function is unbounded, which implies that the beam would carry infinite energy. 18 This effect becomes more severe as the pulse length decreases. On the other hand, in the case of longer (such as 100-fs-long) pulses, Wang's equation approaches the continuous wave Gaussian beam. As the function which describes the continuous Gaussian beam is computationally simpler, we finally decided to use this latter one in the simulator.
Since the publication of our conference paper, 25 we have discovered an error in the former program code: a multiplier factor was missing in the calculation of the illumination intensity from the laser power. This is why we obtained approximately the same excitation probabilities there for laser powers of 30 and 50 mW as here for 25 and 40 mW.

Conclusions
We presented simulation examples which demonstrate quantitatively how the illumination laser power, diffusion, and the photodynamic parameters of the fluorophore (absorption cross-section of photobleaching transitions, lifetime of intersystem crossing, and transition T 1 → S 0 ) affect the number of fluorescent photons and the three-dimensional distribution of fluorescence. We demonstrated how photobleaching changes the spatial distribution of photon emission, producing a "dark hole" in the focus: the simulation results revealed that as more and more laser pulses hit the sample and photobleaching progresses, the most photons are not emitted from the focus but from an ellipsoid situated around the focus. Increased laser power, higher photobleaching cross-section, as well as longer illumination can produce this "dark hole." On the other hand, diffusion can inhibit this phenomenon. Our results evidenced that photobleaching not only decreases the number of emitted photons, but also deteriorates the resolution of the fluorescent techniques.