**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

**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

#### Open Access

## Abstract

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}^{–}^{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}^{–}^{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.

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.

**F1 :**

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$.

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):

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

^{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

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

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}

In the following, we recall the theoretical background of one- and 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:

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

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 $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:

- $S0\u2192S1$, $S1\u2192B2$, and $T1\u2192B4$ via two-photon absorption characterized by absorption cross-sections $\delta exc$, $\delta pb2$, and $\delta pb4$,
- $S1\u2192B1$ and $T1\u2192B3$ via one-photon absorption characterized by absorption cross-sections $\sigma pb1$ and $\sigma pb3$,
- $S1\u2192S0$ radiative relaxation followed by emission of a fluorescent photon, characterized by time constant $\tau f$,
- $S1\u2192S0$ nonradiative relaxation (internal conversion) characterized by time constant $\tau IC$,
- $S1\u2192T1$ intersystem crossing characterized by time constant $\tau ISC$, and
- $T1\u2192S0$ nonradiative relaxation characterized by time constant $\tau TS$.

**F2 :**

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 $A\u2192B$) is

- 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- ($\sigma pb1$ and $\sigma pb3$) or two-photon absorption cross-section ($\delta exc$, $\delta pb2$, and $\delta pb4$) which characterizes the process. The rate constant can be calculated as
Display Formula $kAB(r,z,t)=\sigma AB\lambda exchcI(r,z,t),$(21)if the transition is induced by one-photon absorption, and asDisplay Formula $kAB(r,z,t)=12\delta AB(\lambda exchc)2I2(r,z,t),$(22)if it is coupled to two-photon absorption,^{11}where $\sigma AB$ and $\delta AB$ are the one- and two-photon absorption cross-sections at wavelength $\lambda 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. - 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/\tau f$, $kIC=1/\tau IC$, $kISC=1/\tau ISC$, and $kTS=1/\tau TS$).

We suppose that the pulse duration is much shorter than the time constant of the relaxation processes (i.e., $\tau pulse\u226a\tau IC$, $\tau f$, $\tau ISC$, $\tau TS$), furthermore the pulse length is much shorter than the pulse repetition time (i.e., $\tau pulse\u226a\tau 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}^{–}^{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:

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.

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 $\u22124\sigma pulse$ to $+4\sigma pulse$, where $\sigma 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 $\tau step=\tau 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=\tau 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 $S1\u2192T1\u2192S0$ 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.

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.

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

^{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

^{23}Separating the diffusion along the three axes, we obtain

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 $\tau step=\tau 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=\u2212\Delta 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=\Delta x$: the dislocation of the molecule is one volume cell toward the positive direction along the axis.

- $p1=Pr(X=x1)=Pr(X=\u2212\Delta x)=p$,
- $p2=Pr(X=x2)=Pr(X=0)=1\u22122p$, and
- $p3=Pr(X=x3)=Pr(X=\Delta x)=p$,

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

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\u2192\u221e$, $Y(k)$ is asymptotically normal (Gaussian) with mean

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

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 $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:

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.

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

^{19}In the case of full aperture detection,

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,

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.

**F4 :**

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 $\tau pulse$. The time period between two consecutive pulses is $\tau rep$. Note that in the figure the ratio of these two periods is distorted for the sake of demonstrativeness; in fact, $\tau pulse$ is much shorter than $\tau 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.

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 $\tau pulse=100\u2009\u2009fs$, pulse repetition time $\tau rep=12.5\u2009\u2009ns$ (i.e., the pulse repetition rate was 80 MHz), and the wavelength of the illuminating light $\lambda exc=800\u2009\u2009nm$, 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 $\tau step=\tau rep$. We set $\theta =60\u2009\u2009deg$ for the beam divergence. In every simulation, we used the arbitrarily chosen $1\u2009\u2009\mu 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\xd710\u221249\u2009\u2009cm4\u2009s=20\u2009\u2009GM$ for the excitation two-photon absorption cross-section $\delta exc$ and 3 ns for the lifetime $\tau 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 $\tau IC=\u221e$ and $\tau ISC=\u221e$. (The one exception is indicated later, see Sec. 3.2.5.)

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 $\sigma pb1$, $\delta pb2$, $\sigma pb3$, and $\delta 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\u2009\u2009\mu m$ and $H=10\u2009\u2009\mu m$, respectively, and the cell size to $\Delta x=0.01\u2009\u2009\mu m$. We used $T=1.25\xd710\u22126\u2009\u2009s$ 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.

**F5 :**

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).

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\xd710\u22126\u2009\u2009s$ on a small region ($R=0.32\u2009\u2009\mu m$ and $H=0.96\u2009\u2009\mu 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 $\Delta x=0.01\u2009\u2009\mu 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 $\Delta x=1.25\xd710\u22123\u2009\u2009\mu m$ as the reference.]

Here, we investigate how much it modifies the emission profile if the time step $\tau 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=25\u2009\u2009mW$ and $\sigma pb1=5\xd710\u221222\u2009\u2009cm2$ 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\xd710\u221210\u2009\u2009m2/s$. The simulation length was $T=1.25\xd710\u22123\u2009\u2009s$ (100,000 cycles), and we set the radius and the height of the simulation volume to $R=0.5\u2009\u2009\mu m$ and $H=1.5\u2009\u2009\mu 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,\u2026,128$. In other words, the simulation time step $\tau step$ was set to the pulse repetition time $\tau rep=12.5\u2009\u2009ns$ of the mode-locked laser and its $1/2,1/4,\u2026,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 $\tau step$ tends to zero.

**F8 :**

(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\u2009\u2009\mu m$ and $H=1.5\u2009\u2009\mu 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\u2009\u2009\mu m$ to $3.0\u2009\u2009\mu 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\u2009\u2009\mu m$ and height $H=1.5\u2009\u2009\mu 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\u2009\u2009\mu m$ and height $H=9\u2009\u2009\mu 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\u2009\u2009\mu m$ and $4.5\u2009\u2009\mu 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\u2009\u2009\mu m$ and $4.5\u2009\u2009\mu m$) in all those simulation cases when diffusion is enabled.

**F9 :**

(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\u2009\u2009\mu m$, $1.0\u2009\u2009\mu m$, and $1.5\u2009\u2009\mu 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\u2009\u2009\mu m$ and height $H=1.5\u2009\u2009\mu 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\u2009\u2009\mu m$ and height $H=9\u2009\u2009\mu m$). Note the difference between the scales of (a)–(c).

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\xd710\u221210\u2009\u2009m2/s$), and for photobleaching cross-section $\sigma pb1$: $5\xd710\u221223\u2009\u2009cm2$ (a value found in the paper of Niesner et al.^{10}) and a value 10 times larger ($5\xd710\u221222\u2009\u2009cm2$) 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.

**F10 :**

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 $\sigma 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\u2009\u2009\mu m$ and $H=4.5\u2009\u2009\mu m$, respectively; however, only the focal part of the profile is shown here.

Even when both the laser power and photobleaching cross-section $\sigma pb1$ were moderate (25 mW and $5\xd710\u221223\u2009\u2009cm2$, middle row of the table), the decline of fluorescence due to photobleaching was perceptible after a $1.25\xd710\u22124-s$-long illumination [10,000 laser cycles, Fig. 10(f)], and it became severe after $1.25\xd710\u22123\u2009\u2009s$ [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.)

###### 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 $\sigma pb1=5\xd710\u221223\u2009\u2009cm2$ again, whereas we set $\delta pb2$ to $6.534\xd710\u221254\u2009\u2009cm4\u2009s$ 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=40\u2009\u2009mW$ 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.

**F12 :**

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\u2009\u2009\mu 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\u2009\u2009\mu 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.

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 $\sigma pb1=5\xd710\u221223\u2009\u2009cm2$ and $\delta pb2=6.534\xd710\u221254\u2009\u2009cm4\u2009s$ as references, and we swept these parameters in the range of three orders of magnitude. The results are summarized in Fig. 13.

**F13 :**

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) $\sigma pb1=5\xd710\u221223\u2009\u2009cm2$ and (b) $\delta pb2=6.534\xd710\u221254\u2009\u2009cm4\u2009s$ (denoted by larger squares). (The radius and height of the simulated region were $R=1.5\u2009\u2009\mu m$ and $H=4.5\u2009\u2009\mu m$, respectively.)

When we decreased $\sigma pb1$ by a factor of 32, the number of emitted photons increased by only 28%, whereas when we increased $\sigma 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 $\u221238%$ 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 ($\sigma pb1=5\xd710\u221223\u2009\u2009cm2$ and $\delta pb2=6.534\xd710\u221254\u2009\u2009cm4\u2009s$); the reference photobleaching cross-section values were both doubled (darker curves) and halved (lighter curves) five times.

**F14 :**

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) $\sigma pb1=5\xd710\u221223\u2009\u2009cm2$ and (b) $\delta pb2=6.534\xd710\u221254\u2009\u2009cm4\u2009s$. 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\u2009\u2009\mu m$ and $H=4.5\u2009\u2009\mu m$, respectively.)

To investigate the effect of diffusion on the number of fluorescent photons, we swept the diffusion constant $D$ from 0 to $8\xd710\u221210\u2009\u2009m2/s$. (For the sake of comparison, the diffusion constant of glycine in water is about $10\u22129\u2009\u2009m2/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 ($\sigma pb1=5\xd710\u221223\u2009\u2009cm2$), 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.

In the simulations presented up to this point, we disabled intersystem crossing by setting its lifetime to $\u221e$. 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 ($\tau ISC$) is, the fewer photons were emitted, especially if transition $T1\u2192S0$ was slow, i.e., $\tau 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.

**F17 :**

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,\u2026,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\u2009\u2009\mu m$ and $H=4.5\u2009\u2009\mu m$, respectively.)