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

Imre B. Juhász; Árpád I. I. Csurgay

doi:10.1117/1.JBO.20.1.015001

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

Imre B. Juhász, Árpád I. I. Csurgay

Author Affiliations +

Journal of Biomedical Optics, Vol. 20, Issue 1, 015001 (January 2015). https://doi.org/10.1117/1.JBO.20.1.015001

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

1. Introduction

Fluorescence is a two-phase process during which photon absorption is followed by photon emission. Its high sensitivity and specificity¹ 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)²^,³ and fluorescence recovery after photobleaching (FRAP) (also known as fluorescence photobleaching recovery) measurements.⁴ 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.⁵ 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⁶) 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).⁷^,⁸

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,⁹ and it also causes considerable difficulties in FCS;²^,³ 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,⁹^–¹¹ and some of them also lack the consideration of photobleaching.¹² 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.²^,³^,¹³^–¹⁵

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.

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

Fig. 1

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

2.1.

Model of the Illumination

A monochromatic, continuous wave laser beam is commonly described as a Gaussian beam (also called Gaussian–Lorentzian profile,⁹ which is a solution of the paraxial Helmholtz equation).¹⁶ 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,⁶ 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,¹⁷ 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):

Eq. (1)

P (t) = P_{\max} \exp (- \frac{t^{2}}{2 σ_{pulse}^{2}}) .

In which

Eq. (2)

σ_{pulse} = \frac{τ_{pulse}}{2 \sqrt{2 \ln 2}}

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

Eq. (3)

P_{\max} = P_{avg} 2 \sqrt{\frac{\ln 2}{π}} \frac{τ_{rep}}{τ_{pulse}},

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

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

Eq. (4)

I (r, z) = I_{0} {[\frac{W_{E 0}}{W_{E} (z)}]}^{2} \exp [- \frac{2 r^{2}}{W_{E}^{2} (z)}],

where

Eq. (5)

W_{E} (z) = W_{E 0} \sqrt{1 + {(\frac{z}{z_{E 0}})}^{2}},

in which

Eq. (6)

z_{E 0} = \frac{π W_{E 0}^{2}}{λ_{exc}}

is the Rayleigh range,

Eq. (7)

W_{E 0} = \frac{λ_{exc}}{π θ},

is the waist radius,

λ_{exc}

is the wavelength of the excitation light, and

θ

is the beam divergence.¹⁶ 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

Eq. (8)

I_{0} = \frac{2 P}{π W_{E 0}^{2}} .

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

Eq. (9)

I (r, z, t) = I_{0} (t) {[\frac{W_{E 0}}{W_{E} (z)}]}^{2} \exp [- \frac{2 r^{2}}{W_{E}^{2} (z)}],

moreover

Eq. (10)

I_{0} (t) = \frac{2 P (t)}{π W_{E 0}^{2}} .

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.¹⁹

2.2.

Photodynamic Model of the Fluorophore Molecules

2.2.1.

Theoretical description of the rates of one- and two-photon absorptions

In the following, we recall the theoretical background of one- and two-photon absorptions according to the books of Boyd²⁰ and Masters and So¹ 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:

Eq. (11)

[\frac{{\vec{p}}^{2}}{2 m} + V (\vec{r})] Ψ (\vec{r}, t) = i ℏ \frac{\partial Ψ (\vec{r}, t)}{\partial t},

whose eigenfunctions and eigenvalues give the possible states and energy levels of the system, where

Ψ

is the state function,

\vec{p}

is the momentum,

m

is the mass of the particle,

V (\vec{r})

is the potential energy at point

\vec{r}

,

i

is the imaginary unit, and

ℏ

is the reduced Planck constant. According to the electric dipole approximation, illumination of the system can be interpreted as a perturbation during which the electric field induces electric dipoles. Consequently, a new term is to be added to the Schrödinger equation when the atomic system is illuminated:

Eq. (12)

[\frac{{\vec{p}}^{2}}{2 m} + V (\vec{r}) - \sum e_{j} \vec{E} (ω, t) \cdot {\vec{r}}_{j}] Ψ (\vec{r}, t) = i ℏ \frac{\partial Ψ (\vec{r}, t)}{\partial t},

where

{\vec{r}}_{j}

is a vector connecting the separated charge

e_{j}

in the

j

’th electric dipole induced by electric field

\vec{E} (ω, t)

of angular frequency

ω

at time

t

. Assuming that the illuminating light is monochromatic and linearly polarized, from Eq. (12) applying the first-order time-dependent perturbation theory, one can express the rate of the transition from the ground state

m

to an excited state

n

due to one-photon absorption as

Eq. (13)

R_{m n}^{(1)} = σ_{m n}^{(1)} (ω) I,

where

σ_{m n}^{(1)} (ω)

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 second-order time-dependent perturbation theory results in the following rate for the transition from the ground state

m

to an excited state

n

, via a virtual state

k

:

Eq. (14)

R_{m n}^{(2)} = σ_{m n}^{(2)} (ω) I^{2},

where

σ_{m n}^{(2)} (ω)

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:

Eq. (15)

σ_{m n}^{(1)} (ω) = \frac{4 π^{2} ω}{n c} | {\vec{μ}}_{m n} | ρ (ω = ω_{m n})

and

Eq. (16)

σ_{m n}^{(2)} (ω) = \frac{8 π^{3} ω^{2}}{n^{2} c^{2} ℏ} {| \sum_{k} \frac{{\vec{μ}}_{m k} {\vec{μ}}_{k n}}{ω_{k n} - ω} |}^{2} ρ (2 ω = ω_{m n}),

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

{\vec{μ}}_{m n}

, which can also be calculated theoretically, namely from the wave functions

φ_{m}

and

φ_{n}

of the two states:

Eq. (17)

{\vec{μ}}_{m n} = \sum_{j} e_{j} \int φ_{m}^{* (0)} {\vec{r}}_{j} φ_{n}^{(0)} d^{3} r .

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

Eq. (18)

I = \frac{n c}{2 π ℏ ω} {| \vec{E} |}^{2} .

2.2.2.

Photodynamic states and transitions in our model

We model the states and state transitions by a multilevel photodynamic model taken from the literature⁹^–¹¹ 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:

• $S_{0} \to S_{1}$ , $S_{1} \to B_{2}$ , and $T_{1} \to B_{4}$ via two-photon absorption characterized by absorption cross-sections $δ_{exc}$ , $δ_{p b 2}$ , and $δ_{p b 4}$ ,
• $S_{1} \to B_{1}$ and $T_{1} \to B_{3}$ via one-photon absorption characterized by absorption cross-sections $σ_{p b 1}$ and $σ_{p b 3}$ ,
• $S_{1} \to S_{0}$ radiative relaxation followed by emission of a fluorescent photon, characterized by time constant $τ_{f}$ ,
• $S_{1} \to S_{0}$ nonradiative relaxation (internal conversion) characterized by time constant $τ_{IC}$ ,
• $S_{1} \to T_{1}$ intersystem crossing characterized by time constant $τ_{ISC}$ , and
• $T_{1} \to S_{0}$ nonradiative relaxation characterized by time constant $τ_{TS}$ .

Fig. 2

The Jablonski diagram of the fluorophore represents the photodynamic states and the possible state transitions of the model. $S_{0}$ , $S_{1}$ , and $T_{1}$ denote the ground state, excited state, and triplet state, respectively; $B_{1}$ , $B_{2}$ , $B_{3}$ , and $B_{4}$ 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.⁹^,¹¹^,¹⁹ 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 \to B$ ) is

Eq. (19)

\dot{A} (r, z, t) = - k_{A B} (r, z, t) A (r, z, t),

whereas the time derivative of the number of molecules in state

B

(supposing the change is due only to transition

A \to B

) is

Eq. (20)

\dot{B} (r, z, t) = k_{A B} (r, z, t) A (r, z, t),

where

k_{A B}

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_{p b 1}$ and $k_{p b 3}$ ) or quadratically ( $k_{exc}$ , $k_{p b 2}$ , and $k_{p b 4}$ ), and it is accordingly proportional to either the one- ( $σ_{p b 1}$ and $σ_{p b 3}$ ) or two-photon absorption cross-section ( $δ_{exc}$ , $δ_{p b 2}$ , and $δ_{p b 4}$ ) which characterizes the process. The rate constant can be calculated as
Eq. (21)
$k_{A B} (r, z, t) = σ_{A B} \frac{λ_{exc}}{h c} I (r, z, t),$
if the transition is induced by one-photon absorption, and as
Eq. (22)
$k_{A B} (r, z, t) = \frac{1}{2} δ_{A B} {(\frac{λ_{exc}}{h c})}^{2} I^{2} (r, z, t),$
if it is coupled to two-photon absorption,¹¹ where $σ_{A B}$ and $δ_{A B}$ are the one- and two-photon absorption cross-sections at wavelength $λ_{exc}$ , respectively, $h$ is Planck’s constant, $c$ is the velocity of light, whereas $I (r, z, t)$ is the illuminating light intensity at point $(r, z)$ at time instance $t$ in units of $W / {cm}^{2}$ . The space- and time-dependences of the rate constants arise from the space- and time-dependence of the illumination intensity. In Eq. (22), the factor $1 / 2$ expresses that two photons are required for the two-photon excitation process.
2. The rest of the transitions (relaxation transitions) are intensity-independent and thus space- and time-independent. Each of them is characterized by a global time constant. The rate constant of such a process is determined as the inverse of the time constant of the process ( $k_{f} = 1 / τ_{f}$ , $k_{IC} = 1 / τ_{IC}$ , $k_{ISC} = 1 / τ_{ISC}$ , and $k_{TS} = 1 / τ_{TS}$ ).

We suppose that the pulse duration is much shorter than the time constant of the relaxation processes (i.e., $τ_{pulse} ≪ τ_{IC}$ , $τ_{f}$ , $τ_{ISC}$ , $τ_{TS}$ ), furthermore the pulse length is much shorter than the pulse repetition time (i.e., $τ_{pulse} ≪ τ_{rep}$ ). Therefore, during a light pulse, the relaxation processes can be neglected and it can be assumed that a fluorophore molecule is excited atmost once per pulse.⁹^–¹¹^,¹⁹ 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).¹⁰ Thus, the following two differential equation systems have to be constructed for each volume element:

Eq. (23)

[\begin{matrix} {\dot{S}}_{0} \\ {\dot{S}}_{1} \\ {\dot{T}}_{1} \\ {\dot{B}}_{1} \\ {\dot{B}}_{2} \\ {\dot{B}}_{3} \\ {\dot{B}}_{4} \end{matrix}] = [\begin{matrix} - k_{exc} & 0 & 0 \\ k_{exc} & - (k_{p b 1} + k_{p b 2}) & 0 \\ 0 & 0 & - (k_{p b 3} + k_{p b 4}) \\ 0 & k_{p b 1} & 0 \\ 0 & k_{p b 2} & 0 \\ 0 & 0 & k_{p b 3} \\ 0 & 0 & k_{p b 4} \end{matrix}] \cdot [\begin{matrix} S_{0} \\ S_{1} \\ T_{1} \end{matrix}]

and

Eq. (24)

[\begin{matrix} {\dot{S}}_{0} \\ {\dot{S}}_{1} \\ {\dot{T}}_{1} \end{matrix}] = [\begin{matrix} (k_{IC} + k_{f}) & k_{TS} \\ - (k_{IC} + k_{f} + k_{ISC}) & 0 \\ k_{ISC} & - k_{TS} \end{matrix}] \cdot [\begin{matrix} S_{1} \\ T_{1} \end{matrix}] .

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.

2.2.3.

Computation of the transition probabilities

Based on the ordinary differential equation system formulated in Eq. (23), one can calculate the probability of each photon-absorption-induced state transition during one laser pulse for every volume cell. Because the pulse train consists of identical laser pulses, it is enough to calculate these probabilities once. $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} \to T_{1} \to S_{0}$ is a two-step decay, for which Bateman gave a solution.²¹]

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

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.

2.3.

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.

2.3.1.

Theoretical description of diffusion: Fick’s law

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

Eq. (25)

J = - D grad ρ,

where

J

is the vector of particle flux density,

D

is the diffusion constant, and

ρ

is the concentration of the particles.²³ The relation between the diffusion flux and the concentration is given by the continuity equation, which expresses the conservation of the particles:

Eq. (26)

div J = - \frac{\partial ρ}{\partial t} .

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

Eq. (27)

\frac{\partial ρ}{\partial t} = D Δ ρ .

The solution of Eq. (27) is

Eq. (28)

ρ (x, y, z, t) = \frac{1}{{(\sqrt{4 π D t})}^{3}} \exp (- \frac{x^{2} + y^{2} + z^{2}}{4 D t}),

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

.²³ Separating the diffusion along the three axes, we obtain

Eq. (29)

ρ (x, y, z, t) = \frac{1}{\sqrt{4 π D t}} \exp (- \frac{x^{2}}{4 D t}) \cdot \frac{1}{\sqrt{4 π D t}} \exp (- \frac{y^{2}}{4 D t}) \cdot \frac{1}{\sqrt{4 π D t}} \exp (- \frac{z^{2}}{4 D t}),

i.e., the concentration of the molecules along each axis is described by a Gaussian distribution with zero mean and variance of

2 D t

.

2.3.2.

Discretized diffusion model

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

• $x_{1} = - Δ x$ : the dislocation of the molecule is one volume cell toward the negative direction along the axis,
• $x_{2} = 0$ : the molecule stays at the same volume cell, and
• $x_{3} = Δ x$ : the dislocation of the molecule is one volume cell toward the positive direction along the axis.

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

• $p_{1} = \Pr (X = x_{1}) = \Pr (X = - Δ x) = p$ ,
• $p_{2} = \Pr (X = x_{2}) = \Pr (X = 0) = 1 - 2 p$ , and
• $p_{3} = \Pr (X = x_{3}) = \Pr (X = Δ x) = p$ ,

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

Eq. (30)

E [X] = \sum_{i = 1}^{3} p_{i} x_{i} = 0

and

Eq. (31)

Var [X] = E [{(X - E [X])}^{2}] = \sum_{i = 1}^{3} p_{i} x_{i}^{2} = 2 p \cdot {(Δ x)}^{2} .

Let $Y^{(k)}$ be a random variable which denotes the dislocation of a molecule after $k$ time steps, i.e.,

Eq. (32)

Y^{(k)} = \binom{\underset{︸}{X + X + \dots + X}}{k} = k \cdot X .

According to the central limit theorem,¹⁷ as $k \to \infty$ , $Y^{(k)}$ is asymptotically normal (Gaussian) with mean

Eq. (33)

E [Y^{(k)}] = μ = \sum_{i = 1}^{k} E [X] = k \cdot E [X] = 0

and variance

Eq. (34)

Var [Y^{(k)}] = σ^{2} = \sum_{i = 1}^{k} Var [X] = k \cdot Var [X] = 2 k \cdot p \cdot {(Δ x)}^{2} .

Equation (29) states that the distribution of the diffusing matter is Gaussian with zero mean and a variance of $2 D t$ ; thus, for the variance $σ^{2}$ of $Y^{(k)}$ , equation

Eq. (35)

σ^{2} = 2 k \cdot p \cdot {(Δ x)}^{2} = 2 D t,

should hold. Substituting

t = k τ_{step}

into Eq. (35), we obtain

Eq. (36)

p = \frac{D τ_{step}}{Δ x^{2}},

for the probability that a molecule moves to one of the neighboring volume cells along a coordinate axis during one time step.

2.3.3.

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.

Fig. 3

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

Because of the cylindrical symmetry, the concentration at point $P$ has to be equal to the concentration at point $Q$ since both points are at distance $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:

Eq. (37)

c_{P} = c_{Q} \approx c_{A} + (r_{0} - ⌊ r_{0} ⌋) (c_{B} - c_{A}),

where

⌊ \dots ⌋

denotes the floor function. Because the distance between points

A

and

B

is unity,

⌊ r_{0} ⌋

is equal to the distance between points

O

and

A

.

2.3.4.

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.

2.4.

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

Eq. (38)

p_{\det} (r, z) = η \cdot Φ \cdot Y (r, z),

where

η

is the efficiency of the detection,

Φ

is the fractional solid angle of the observation, and

Y (r, z)

is the observation beam profile.¹⁹ In the case of full aperture detection,

Eq. (39)

Φ = \frac{{NA}^{2}}{4 n^{2}}

and

Eq. (40)

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,

Eq. (41)

Φ = \frac{λ_{emit}^{2}}{4 π^{2} n^{2} W_{D 0}^{2}}

and

Eq. (42)

Y (r, z) = {[\frac{W_{D 0}}{W_{D} (z)}]}^{2} \exp [- \frac{2 r^{2}}{W_{D}^{2} (z)}];

i.e., the observation beam profile is described by a Gaussian beam, whose waist radius is

W_{D 0}

.

λ_{emit}

denotes the wavelength of the emitted fluorescent light which is to be detected. (Subscript

D

refers to detection.).

2.5.

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.

Fig. 4

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.

3. Simulations and Results

3.1.

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¹⁰^,¹¹ for real fluorophores. We used $2 \times 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} = \infty$ and $τ_{ISC} = \infty$ . (The one exception is indicated later, see Sec. 3.2.5.)

3.1.1.

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 $σ_{p b 1}$ , $δ_{p b 2}$ , $σ_{p b 3}$ , and $δ_{p b 4}$ 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 \times 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.

Fig. 5

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

Fig. 6

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

3.1.2.

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 \times 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 = 0.01 μ m$ for the cell size in further simulations. [For this cell size, the error of the emitted-photon number is 0.60% and 0.69% for the two chosen laser power values (25 and 40 mW) considering the case with $Δ x = 1.25 \times 10^{- 3} μ m$ as the reference.]

Fig. 7

Dependence of (a) the fluorescent photon number and (b) the number of simulated cells on the cell size. In (a), the reference is the photon number obtained for the simulation case with $Δ x = 1.25 \times 10^{- 3} μ m$ .

3.1.3.

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 $σ_{p b 1} = 5 \times 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 \times 10^{- 10} m^{2} / s$ . The simulation length was $T = 1.25 \times 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, \dots, 128$ . In other words, the simulation time step $τ_{step}$ was set to the pulse repetition time $τ_{rep} = 12.5 ns$ of the mode-locked laser and its $1 / 2, 1 / 4, \dots, 1 / 128$ .

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

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

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

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

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

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

3.2.

Investigation of the Three-Dimensional Fluorescent Profile and Its Time Evolution

3.2.1.

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 \times 10^{- 10} m^{2} / s$ ), and for photobleaching cross-section $σ_{p b 1}$ : $5 \times 10^{- 23} {cm}^{2}$ (a value found in the paper of Niesner et al.¹⁰) and a value 10 times larger ( $5 \times 10^{- 22} {cm}^{2}$ ) 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.

Fig. 10

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 $σ_{p b 1}$ were moderate. Even in this case, the decline of fluorescence due to photobleaching was perceptible after an illumination of 10,000 laser cycles (f), and it became severe after 100,000 cycles (g), especially at the focus, where a “dark hole” evolved. Laser power was increased in the first row [(a)–(d)] and the photobleaching cross-section in the third one [(i)–(l)]. In both cases, photobleaching was accelerated. In columns 1–3, diffusion was disabled; however, it occurred in column 4 compensating photobleaching more or less. Note that in the simulations that are presented in column 4 and in the case for which diffusion was enabled, the radius and height of the simulated region were $R = 1.5 μ m$ and $H = 4.5 μ m$ , respectively; however, only the focal part of the profile is shown here.

Even when both the laser power and photobleaching cross-section $σ_{p b 1}$ were moderate (25 mW and $5 \times 10^{- 23} {cm}^{2}$ , middle row of the table), the decline of fluorescence due to photobleaching was perceptible after a $1.25 \times 10^{- 4} - s$ -long illumination [10,000 laser cycles, Fig. 10(f)], and it became severe after $1.25 \times 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.

Fig. 11

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

3.2.2.

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 $σ_{p b 1} = 5 \times 10^{- 23} {cm}^{2}$ again, whereas we set $δ_{p b 2}$ to $6.534 \times 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 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.

Fig. 12

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

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.

3.2.3.

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 $σ_{p b 1} = 5 \times 10^{- 23} {cm}^{2}$ and $δ_{p b 2} = 6.534 \times 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.

Fig. 13

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) $σ_{p b 1} = 5 \times 10^{- 23} {cm}^{2}$ and (b) $δ_{p b 2} = 6.534 \times 10^{- 54} {cm}^{4} s$ (denoted by larger squares). (The radius and height of the simulated region were $R = 1.5 μ m$ and $H = 4.5 μ m$ , respectively.)

When we decreased $σ_{p b 1}$ by a factor of 32, the number of emitted photons increased by only 28%, whereas when we increased $σ_{p b 1}$ by the same factor, the number of emitted photons fell to 50% of the reference.

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

Figure 14 shows the decay of photon emission in time due to photobleaching. The thick curves belong to the aforementioned reference cases ( $σ_{p b 1} = 5 \times 10^{- 23} {cm}^{2}$ and $δ_{p b 2} = 6.534 \times 10^{- 54} {cm}^{4} s$ ); the reference photobleaching cross-section values were both doubled (darker curves) and halved (lighter curves) five times.

Fig. 14

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) $σ_{p b 1} = 5 \times 10^{- 23} {cm}^{2}$ and (b) $δ_{p b 2} = 6.534 \times 10^{- 54} {cm}^{4} s$ . These values were both doubled and halved five times: the higher the photobleaching cross-section is, the lower and darker the curve is. (The radius and height of the simulated region were $R = 1.5 μ m$ and $H = 4.5 μ m$ , respectively.)

3.2.4.

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 \times 10^{- 10} m^{2} / s$ . (For the sake of comparison, the diffusion constant of glycine in water is about $10^{- 9} m^{2} / s$ .²⁴) 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 ( $σ_{p b 1} = 5 \times 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.

Fig. 15

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

3.2.5.

Effect of the intersystem crossing

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

The simulation results (Fig. 16) corresponded to the expectations. The smaller the lifetime of intersystem crossing ( $τ_{ISC}$ ) is, the fewer photons were emitted, especially if transition $T_{1} \to S_{0}$ was slow, i.e., $τ_{TS}$ was large. The return from the triplet state $T_{1}$ to the ground state $S_{0}$ 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.

Fig. 16

Dependence of the number of emitted photons (a) on the lifetime of intersystem crossing ( $τ_{ISC}$ ) and (b) on the lifetime of the transition $T_{1} \to S_{0}$ ( $τ_{TS}$ ). (The radius and height of the simulated region were $R = 1.5 μ m$ and $H = 4.5 μ m$ , respectively.)

Fig. 17

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 $1 st, 50 th, 100 th, \dots, 2000 th$ 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 $T_{1}$ over time. (The radius and height of the simulated region were $R = 1.5 μ m$ and $H = 4.5 μ m$ , respectively.)

If intersystem crossing is permitted, then the fluorophore molecules in the triplet state $T_{1}$ might be photobleached by the absorption of one or two additional photons. According to the simulations of this scenario (not presented in detail), these photobleaching transitions affect the photon emission profile similar to those that start from the excited singlet state $S_{1}$ (see Figs. 10 to 12). Nevertheless, in the former case, the evolution of the photon distribution is also driven by intersystem crossing, which flattens the profile with the “dark hole” at the focus, formed by photobleaching.

4. Discussion

Although our model is founded on the physical principles presented in Sec. 2, its validation and calibration require experimental measurements. Besides the approximations and the neglected factors indicated in the description of the model, the most remarkable difference between the presented simulations and real 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-per-cycle 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 two-photon 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.²⁵ At that stage of our work, we used an equation given by Wang et al.²⁶ 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.¹⁸ 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,²⁵ 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.

5. 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} \to 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.

Acknowledgments

We, the authors, gratefully acknowledge the grant from TÁMOP-4.2.1.B-11/2/KMR-2011-0002. We are grateful to the Multidisciplinary Doctoral School of Sciences and Technology, Pázmány Péter Catholic University, Budapest, Hungary and its leaders, Professor Tamás Roska and Professor Péter Szolgay for providing the conditions of the research work. We thank also our colleague, Ádám Fekete for the inspiring discussions.

Appendices

Appendix:

Pseudocode of the Algorithm

% calculation of probabilities of transitions and diffusion during one simulation step
for each volume cell $(r, z)$ do
for each photon-absorption-induced state transition $A \to B$ do
$P r_{abs, A \to B} (r, z)$ ← numerical solution of Eq. (23) by MATLAB function ode45
end for
for each relaxation transition $A \to B$ do
$P r_{relax, A \to B} (r, z)$ ← numerical solution of Eq. (24) by MATLAB function ode45
end for
end for
$P r_{diff}$ ← value given by Eq. (36)
% initial state: every fluorophore molecule is in ground state
for each volume cell $(r, z)$ do
$S_{0} (r, z) \leftarrow N_{init}$
for each other state A do
$A (r, z) \leftarrow 0$
end for
end for
for $i = 1$ → number of laser cycles do
% simulation of photon-absorption-induced transitions
for each volume cell $(r, z)$ do
for each photon-absorption-induced state transition $A \to B$ do
$Δ A_{A \to B} (r, z) \leftarrow A \cdot P r_{abs, A \to B} (r, z)$
end for
for each state A do
$A (r, z) \leftarrow A (r, z) + \sum_{\forall X \neq A} [- Δ A_{A \to X} (r, z) + Δ A_{X \to A} (r, z)]$
end for
end for
for $j = 1$ → number of simulation steps per laser cycle do
% simulation of relaxation transitions
for each volume cell $(r, z)$ do
for each relaxation state transition $A \to B$ do
$Δ A_{A \to B} (r, z) \leftarrow A \cdot P r_{relax, A \to B} (r, z)$
end for
for each state A do
$A (r, z) \leftarrow A (r, z) + \sum_{\forall X \neq A} [- Δ A_{A \to X} (r, z) + Δ A_{X \to A} (r, z)]$
end for
end for
% simulation of diffusion along axes $r$ , $z$ , and the third axis
for each volume cell $(r, z)$ do
$A_{stays} (r, z) \leftarrow A \cdot (1 - 2 \cdot P r_{diff})$
$A_{moves} (r, z) \leftarrow A \cdot P r_{diff}$
end for
for each volume cell $(r, z)$ do
$A (r, z) \leftarrow A_{stays} (r, z) + A_{moves} (r - 1, z) + A_{moves} (r + 1, z)$
end for
for each volume cell $(r, z)$ do
$A_{stays} (r, z) \leftarrow A \cdot (1 - 2 \cdot P r_{diff})$
$A_{moves} (r, z) \leftarrow A \cdot P r_{diff}$
end for
for each volume cell $(r, z)$ do
$A (r, z) \leftarrow A_{stays} (r, z) + A_{moves} (r, z - 1) + A_{moves} (r, z + 1)$
end for
for each volume cell $(r, z)$ do
$A_{stays} (r, z) \leftarrow A \cdot (1 - 2 \cdot P r_{diff})$
$A_{moves} (r, z, *)$ ← value given by Eq. (37)
end for
for each volume cell $(r, z)$ do
$A (r, z) \leftarrow A_{stays} (r, z) + 2 \cdot A_{moves} (r, z, *)$
end for
end for
end for

References

1.

B. R. MastersP. T. C. So, “Classical and quantum theory of one-photon and multiphoton fluorescence spectroscopy,” Handbook of Biomedical Nonlinear Optical Microscopy, 91 –152 Oxford University Press, New York (2008). Google Scholar

2.

A. Delonet al., “Photobleaching, mobility, and compartmentalisation: inferences in fluorescence correlation spectroscopy,” J. Fluoresc., 14 (3), 255 –267 (2004). http://dx.doi.org/10.1023/B:JOFL.0000024557.73246.f9 JOFLEN 1053-0509 Google Scholar

3.

Z. PetrášekP. Schwille, “Photobleaching in two-photon scanning fluorescence correlation spectroscopy,” ChemPhysChem, 9 147 –158 (2008). http://dx.doi.org/10.1002/(ISSN)1439-7641 CPCHFT 1439-4235 Google Scholar

4.

P.-C. Cheng, “The contrast formation in optical microscopy,” Handbook of Biological Confocal Microscopy, 162 –206 3rd ed.Springer, New York (2006). Google Scholar

5.

W. DenkJ. H. StricklerW. W. Webb, “Two-photon laser scanning fluorescence microscopy,” Science, 248 73 –76 (1990). http://dx.doi.org/10.1126/science.2321027 SCIEAS 0036-8075 Google Scholar

6.

W. DenkD. W. PistonW. W. Webb, “Multi-photon molecular excitation laser-scanning microscopy,” Handbook of Biological Confocal Microscopy, 535 –549 3rd ed.Springer, New York (2006). Google Scholar

7.

L. SchermellehR. HeintzmannH. Leonhardt, “A guide to super-resolution fluorescence microscopy,” J. Cell Biol., 190 165 –175 (2010). http://dx.doi.org/10.1083/jcb.201002018 JCLBA3 0021-9525 Google Scholar

8.

B. HuangM. BatesX. Zhuang, “Super-resolution fluorescence microscopy,” Annu. Rev. Biochem., 78 993 –1016 (2009). http://dx.doi.org/10.1146/annurev.biochem.77.061906.092014 ARBOAW 0066-4154 Google Scholar

9.

J. Mertz, “Molecular photodynamics involved in multi-photon excitation fluorescence microscopy,” Eur. Phys. J. D, 3 (1), 53 –66 (1998). http://dx.doi.org/10.1007/s100530050148 EPJDF6 1434-6060 Google Scholar

10.

R. NiesnerW. RothK.-H. Gericke, “Photophysical aspects of single-molecule detection by two-photon excitation with consideration of sequential pulsed illumination,” ChemPhysChem, 5 678 –687 (2004). http://dx.doi.org/10.1002/(ISSN)1439-7641 CPCHFT 1439-4235 Google Scholar

11.

C. EggelingA. VolkmerC. A. M. Seidel, “Molecular photobleaching kinetics of rhodamine 6G by one- and two-photon induced confocal fluorescence microscopy,” ChemPhysChem, 6 791 –804 (2005). http://dx.doi.org/10.1002/(ISSN)1439-7641 CPCHFT 1439-4235 Google Scholar

12.

K. Paloet al., “Calculation of photon-count number distributions via master equations,” Biophys. J., 90 2179 –2191 (2006). http://dx.doi.org/10.1529/biophysj.105.066084 BIOJAU 0006-3495 Google Scholar

13.

D. Mazzaet al., “Role of three-dimensional bleach distribution in confocal and two-photon fluorescence recovery after photobleaching experiments,” Appl. Opt., 46 (30), 7401 –7411 (2007). http://dx.doi.org/10.1364/AO.46.007401 APOPAI 0003-6935 Google Scholar

14.

N. Smisdomet al., “Fluorescence recovery after photobleaching on the confocal laser-scanning microscope: generalized model without restriction on the size of the photobleached disk,” J. Biomed. Opt., 16 (4), 046021 (2011). http://dx.doi.org/10.1117/1.3569620 JBOPFO 1083-3668 Google Scholar

15.

E. A. Schnellet al., “Diffusion measured by fluorescence recovery after photobleaching based on multiphoton excitation laser scanning microscopy,” J. Biomed. Opt., 13 (6), 064037 (2008). http://dx.doi.org/10.1117/1.3042274 JBOPFO 1083-3668 Google Scholar

16.

B. E. A. SalehM. C. Teich, Fundamentals of Photonics, 2nd ed.John Wiley & Sons, Hoboken, New Jersey (2007). Google Scholar

17.

G. A. KornT. M. Korn, Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulas for Reference and Review, Dover Publications, Mineola, New York (2000). Google Scholar

18.

A. April, “Ultrashort, strongly focused laser pulses in free space,” Coherence and Ultrashort Pulse Laser Emission, 355 –382 InTech(2010). Google Scholar

19.

M. Kauertet al., “Absolute measurement of molecular two-photon absorption cross-sections using a fluorescence saturation technique,” Opt. Express, 14 8434 –8447 (2006). http://dx.doi.org/10.1364/OE.14.008434 OPEXFF 1094-4087 Google Scholar

20.

R. W. Boyd, Nonlinear Optics, 3rd ed.Academic Press, San Diego (2008). Google Scholar

21.

H. Bateman, “Solution of a system of differential equations occurring in the theory of radioactive transformations,” Proc. Cambridge Philos. Soc., 15 423 –427 (1910). PCPSA4 0008-1981 Google Scholar

22.

“ode45: Solve nonstiff differential equations; medium order method,” (2014) http://www.mathworks.com/help/matlab/ref/ode45.html November ). 2014). Google Scholar

23.

Handbook of Physics, Springer, New York (2002). Google Scholar

24.

R. G. Mortimer, Physical Chemistry, 2nd ed.Academic Press, Burlington, Massachusetts (2000). Google Scholar

25.

I. B. JuhaszA. I. Csurgay, “Simulation-based investigation of the three-dimensional distribution of fluorescence and photobleaching in multi-photon excited samples,” Proc. SPIE, 9131 91310S (2014). http://dx.doi.org/10.1117/12.2052243 PSISDG 0277-786X Google Scholar

26.

Z. Wanget al., “Space-time profiles of an ultrashort pulsed gaussian beam,” IEEE J. Quantum Electron., 33 (4), 566 –573 (1997). http://dx.doi.org/10.1109/3.563385 IEJQA7 0018-9197 Google Scholar

Biography

Imre B. Juhász is a PhD student at the Faculty of Information Technology and Bionics, Pázmány Péter Catholic University, where he received his MSc degree in computer engineering in 2012. His main interest is modeling and simulation of light–matter interaction.

Árpád I. Csurgay received his Dipl-Eng degree in 1959, his PhD degree in 1964 from TU Budapest, and his DSc degree in 1973 from the Hungarian Academy of Sciences. He was affiliated with TU Budapest and University of Notre Dame. Currently, he is a professor of applied physics at the PPCU Budapest. He is a fellow of the IEEE and a member of the Hungarian Academy of Sciences and Academia Europaea.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Imre B. Juhász and Árpád I. I. Csurgay "Fluorescence in two-photon-excited diffusible samples exposed to photobleaching: a simulation-based study," Journal of Biomedical Optics 20(1), 015001 (9 January 2015). https://doi.org/10.1117/1.JBO.20.1.015001

Published: 9 January 2015

Access the abstract

JOURNAL ARTICLE
19 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 1 scholarly publication.

Explore citations on Lens.org

KEYWORDS

Photons

Diffusion

Molecules

Computer simulations

Luminescence

Optical simulations

Absorption

1.

Introduction

2.

Simulation Model

Fig. 1

2.1.

Model of the Illumination

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

2.2.

Photodynamic Model of the Fluorophore Molecules

2.2.1.

Theoretical description of the rates of one- and two-photon absorptions

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

Eq. (16)

Eq. (17)

Eq. (18)

2.2.2.

Photodynamic states and transitions in our model

Fig. 2

Eq. (19)

Eq. (20)

Eq. (21)

Eq. (22)

Eq. (23)

Eq. (24)

2.2.3.

Computation of the transition probabilities

2.3.

Diffusion Model

2.3.1.

Theoretical description of diffusion: Fick’s law

Eq. (25)

Eq. (26)

Eq. (27)

Eq. (28)

Eq. (29)

2.3.2.

Discretized diffusion model

Eq. (30)

Eq. (31)

Eq. (32)

Eq. (33)

Eq. (34)

Eq. (35)

Eq. (36)

2.3.3.

Cylindrical symmetry

Fig. 3

Eq. (37)

2.3.4.

Boundary conditions

2.4.

Model of the Photon Detection

Eq. (38)

Eq. (39)

Eq. (40)

Eq. (41)

Eq. (42)

2.5.

Process of the Simulation

Fig. 4

3.

Simulations and Results

3.1.

Selection of the Simulation Parameters

3.1.1.

Calibration of the laser power