2.1.

Mathematical Task

The goal of the numerical model of aPDT is to determine the local temporal therapy success. The quantity used to describe the therapy success is the survival fraction $B(\overrightarrow{r},t)$ of bacteria in the dental pocket

The receptor $R$ is the part of the bacterial membrane that reacts with the singlet oxygen ${}^1{\mathrm{O}}_2$. It has no biological equivalent and is introduced as a measure for biological activity. A bacteria is considered to be living as long as a certain amount $r_{\min }$ of its receptors $R$ is active. In order to deduce any information of the aPDT success from the receptor concentration $[R]$, the number of active receptors of an individual bacteria is estimated with the binomial distribution. Only those individual bacteria that have at least $r_{\min }$ receptors still active contribute to the survival fraction $B$. The value $B=0$ of the survival fraction $B$ is the criterion to measure the therapy outcome.

Our current modeling task is a two-dimensional (2-D) model of a periodontium phantom model shown in (Fig. 2). aPDT takes place in the area of the dental pocket ${\mathrm{\Omega }}_{\mathrm{aPDT}}$ only. All the substances taking part in aPDT are homogeneously distributed at the beginning of the therapy. In the area of gingiva ${\mathrm{\Omega }}_{\mathrm{ging}}$ and the dentine ${\mathrm{\Omega }}_{\mathrm{dent}}$ only radiative transfer is considered (Fig. 2). The laser radiation is entering the calculation area with a Gaussian intensity profile at the boundaries ${\mathrm{\Gamma }}_{\mathrm{abs},1}$ or ${\mathrm{\Gamma }}_{\mathrm{abs},2}$. Any caustic of the Gaussian beam is neglected.

Fig. 2

Schematic representation of the calculation area. The dimensions of the area and the boundaries for diffusion processes are seen.

In order to determine the receptor concentration $[R](\overrightarrow{r},t)$ in our calculation area, we revert to the Foster et al. and Hu et al. rate equation model.^13,17 As explained there in detail, the concentrations of the three PS states $[S_0]$, $[S_1]$, and $[T]$, and the two of oxygen $[{}^3{\mathrm{O}}_2]$ and $[{}^1{\mathrm{O}}_2]$, change in time due to photochemical processes (absorption, fluorescence, phosphorescence and resonant collision, cytotoxicity, and photobleaching, see Fig. 1). Those processes and the change of oxygen concentration in space due to diffusion processes are expressed mathematically in Eqs. (2)–(6)

In accordance with Liu et al.,²¹ we only include photobleaching of $S_0$ mediated by singlet oxygen ${}^1{\mathrm{O}}_2$ [see Eqs. (2) and (6)].

The calculation area for the rate equation systems (2)–(6) is ${\mathrm{\Omega }}_{\mathrm{aPDT}}$, shown in Fig. 2. The boundary conditions for all concentrations $[X]$ are considered to be isolated, i.e., mass transport through the boundaries is forbidden, which means all the substances remain inside the calculation area ${\mathrm{\Omega }}_{\mathrm{aPDT}}$. For oxygen in the ground state, values at the boundaries with air ${\mathrm{\Gamma }}_{\mathrm{air}}$ are fixed and equal to the saturation concentration ${[{}^3{\mathrm{O}}_2]}_{\mathrm{sat}}$ of oxygen in water at 37°C

The oxygen in the singlet state ${}^1{\mathrm{O}}_2$ reacts with the bacterial receptors $R$ [Eq. (6)] and thereby inactivates them (cytotoxicity in Fig. 1). On the other hand, it has been observed that the bacteria profit from cell repair mechanisms.²⁷ Although the cell repair rate is not significant on the time scale of a successful aPDT treatment, it needs to be included in the model to cover cases of very slow or uncompleted aPDT progresses. In those cases, the therapy outcome is very sensitive to the cell repair rate. Thus, the concentration of receptors $[R]$ regenerates with the rate $U$

The diffusion of bacteria, and thus receptors, is not considered in this model, since bacteria will not be totally destroyed but only devitalized. This does not lead to gradients in bacteria concentration. Diffusion of PS is not considered in Eqs. (2)–(4), because PS is bound to the bacteria membranes. The initial condition for the bacteria receptors is

The local photon density $\rho (\overrightarrow{r},t)$ is required as input in Eqs. (2) and (3) to induce the photochemical reactions and thus reduce the survival fraction $B(\overrightarrow{r},t)$.

The propagation of light, the chemical processes during aPDT, and the diffusion processes inside the treated area are coupled processes as the local radiation intensity enables the chemical reactions and depends on the local optical properties, which are strongly influenced by the concentration of PS. In addition, due to photobleaching, the local PS concentration varies. The photon density $\rho (\overrightarrow{r},t)$ is determined by

2.2.

Radiative Transfer

Due to its fast time scale, the intensity distribution relaxes quickly to its temporary equilibrium and the stationary equation of radiative transfer

This distribution was found to be accurate for biological tissues.²² $\theta$ indicates the scattering angle and $g$ is the anisotropy coefficient. Since the local optical properties change in the volume according to the rate equations, the photon density distribution has to be updated after every change of the PS concentration $[S_0](\overrightarrow{r},t)$. The local optical properties are stored in a grid with side lengths $\mathrm{\Delta }x$ and $\mathrm{\Delta }y$. When a photon crosses the border between two grid cells, the covered path length is subtracted from the original path length and the rest of it is scaled in order to comply with the local optical properties in the current grid cell (Fig. 3). The photon density inside the calculation area is stored and serves as an input parameter for the rate equations at every grid cell.

Fig. 3

Schematic representation of the dental pocket ${\mathrm{\Omega }}_{\mathrm{aPDT}}$. The processes describing interaction between two grid cells are sketched.

2.3.

Antimicrobial Photodynamic Therapy Reaction Kinetics and Diffusion

The system of partial differential Eqs. (2)–(6) and (12) is solved on the same 2-D grid as the Monte Carlo simulations (Fig. 3). In order to reduce calculation time, we decide to choose two different time steps for the reaction kinetics and diffusion terms. The numerical procedure is shown for the reaction diffusion equation of oxygen in its ground state $[{}^3{\mathrm{O}}_2]$, but is analogous for the other concentrations. The total change of oxygen is

The diffusion term $D_{\mathrm{ox}}{\nabla }^2[{}^3{\mathrm{O}}_2]$ is evaluated at the time step $t_n$ using the finite-difference method (Fig. 3). The result is used as a constant diffusion term in the rate equation for oxygen. The rate equation is solved by the MATLAB solver ode23t until the time $t_{n+1}=t_n+\mathrm{\Delta }t$ is reached. At the time $t_{n+1}$, the diffusion term is evaluated as a first prediction. The rate equations are solved again using a linear interpolation between the diffusion terms at $t_n$ and $t_{n+1}$ and more accurate results of the rate equations are received.

3.1.

Homogeneous Model

In order to gain an understanding of how the concentrations in the irradiated area evolve in time, a zero-dimensional model is solved as a first step. The solution of this point system is equivalent to the solution of a system extended in space with homogeneously distributed initial conditions and processing parameters. In such a system, both transport processes and position play no role for the evolution of concentrations. The evolution of the concentrations in time in such a point system is shown for two sets of initial conditions in Figs. 4 and 5. The therapy success is expressed in the survival fraction $B$ of the bacteria. The therapy is successful if the survival fraction $B$ is reduced to zero during the irradiation time, which is the case in Fig. 4. After the receptors $R$ are significantly reduced, some decrease in PS and oxygen concentration can be noticed. This decrease is the result of the oxidation reaction between PS ($S_0$) and oxygen (${}^1{\mathrm{O}}_2$), called photobleaching (see Fig. 1). Photobleaching leads to a reduction of the PS concentration and is followed by higher penetration of laser light into the irradiated area. In the double logarithmic presentation, the comparatively steep decay of receptors $R$ and survival fraction $B$ is very striking. The therapy success is achieved only after a certain irradiation time with a constant intensity. This behavior of a minimum threshold dose for laser radiation has been observed in PDT for cancer treatment experimentally³¹ and we expect to find it in experiments for aPDT as well.

Fig. 4

Solution of the rate Eqs. (2)–(6) and (12) to describe the reaction kinetics for homogeneous distribution of initial values and a photon density constant in time and space. The initial value of oxygen concentration is ${[{}^3{\mathrm{O}}_2]}_0/{[R]}_0=1.2$. The constant photon density is $\rho =2·{10}^6{\mathrm{cm}}^{-3}$. The other parameters are taken from Table 1.

Fig. 5

Solution of the rate Eqs. (2)–(6) and (12) to describe the reaction kinetics for homogeneous distribution of initial values and a photon density constant in time and space. The initial value of oxygen concentration is ${[{}^3{\mathrm{O}}_2]}_0/{[R]}_0=0.8$. The constant photon density is $\rho =2·{10}^6{\mathrm{cm}}^{-3}$. The other parameters are taken from Table 1.

In addition, some significance of the ratio of the initial value of oxygen concentration and receptor concentration ${[{}^3{\mathrm{O}}_2]}_0/{[R]}_0$ can be seen. In Fig. 4, this ratio is ${[{}^3{\mathrm{O}}_2]}_0/{[R]}_0=1.2$ and the therapy success is achieved after 12 s of irradiation. With a too little initial value for oxygen concentration (${[{}^3{\mathrm{O}}_2]}_0/{[R]}_0=0.8$ in Fig. 5), the therapy is not successful ($B>0$) after 100 s. Instead, the reduction of the survival fraction $B$ is limited because the oxygen concentration is reduced to zero before enough bacterial receptors are oxidized. The constant photon density is $\rho =2·{10}^6{\mathrm{cm}}^{-3}$ in both figures. All the other parameters correspond to the parameter seed point, listed in Table 1.

Using scatterplots, we analyzed the therapy outcome $B=0$ after 30 s of irradiation. We found that the initial value of the oxygen concentration ${[{}^3{\mathrm{O}}_2]}_0$ is the most sensitive parameter of the therapy outcome $B=0$ at our seed point (Table 1). Analyzing the homogeneous model around its parameter seed point, we saw that the ratio between theb initial value of oxygen concentration and receptor concentration ${[{}^3{\mathrm{O}}_2]}_0/{[R]}_0$ has to be at least 1 to reduce the receptor concentration significantly. Therefore, we use ${[{}^3{\mathrm{O}}_2]}_0/{[R]}_0=1.2$ as this ratio for all our calculations. The scatterplot for the survival fraction $B(30\mathrm{s})$ after 30 s of irradiation is shown in Fig. 6. We learn from Fig. 6 that successful aPDT is not possible without the presence of oxygen. If enough oxygen is available at the beginning of the treatment, the reduction of the survival fraction $B$ to 0 depends on other parameters. Therefore, aPDT has to be planned carefully with a deep understanding of the underlying physical and chemical processes.

Fig. 6

Sensitivity of the survival fraction $B(30\mathrm{s})$ after 30 s of irradiation on the initial value of oxygen concentration ${[{}^3{\mathrm{O}}_2]}_0/{[R]}_0$. The simulation results are obtained solving the homogeneous model for the aPDT reaction kinetics. The parameters are varied randomly around the seed point in Table 1. The photon density is varied around the value ${\rho }_0=2·{10}^6{\mathrm{cm}}^{-3}$.

3.2.

Therapy Progress

The characteristic result of our model is shown in Fig. 7. The irradiated and thus successfully ($B=0$) treated area is growing in time. The calculation is performed with the parameters given in Table 1. For the survival fraction, no transition is seen between $B=1$ and $B=0$. This is explained by the steep reduction of $B$ already seen in the homogeneous model (Fig. 4). At the positions $\overrightarrow{r}$, where the survival fraction of bacteria $B$ is reduced to virtually zero, photobleaching occurs and photons travel further inside the volume. Therefore, additional to the time required by the reaction kinetics to reduce the survival fraction of bacteria $B$ (12 s in Fig. 4), the time that photons need to reach the position $\overrightarrow{r}$ has to be considered for the therapy. This is visualized in Fig. 8 for the same calculations. Therefore, the evolution of survival fractions $B$ at three different heights of the dental pocket is seen. As expected, the survival fraction is much more quickly reduced on the top of the dental pocket ${\mathrm{\Omega }}_{\mathrm{aPDT}}$ than at its bottom. There, after some seconds only, a sufficient photon density is reached.

Fig. 7

Photon density $\rho (\overrightarrow{r},t)$ and survival fraction $B(\overrightarrow{r},t)$ at five points in time inside the dental pocket ${\mathrm{\Omega }}_{\mathrm{aPDT}}$. The parameters for the calculation are given in Table 1.

Fig. 8

Survival fraction of bacteria $B(\overrightarrow{r},t)$ at three different heights of the dental pocket. The computation parameters are given in Table 1.

Photobleaching and its influence on the optical properties in tissue have been examined experimentally and theoretically for PDT in the treatment of tumors.^18,32 In these studies, fluorescence signals of PS have been detected over time. Their dependence on the applicated laser dose supports our idea of increasing the successfully treated area.

In order to compare different settings of parameters, we introduce the quantities successfully treated area $A_{\mathrm{success}}$

Fig. 9

Ratio of successfully treated and total area $A_{\mathrm{success}}/A_0$ evolving in time for four different sets of initial conditions. The initial conditions different from the parameter seed point (Table 1) are labeled in the graph.

Our simulation shows a faster aPDT success for the PS initial concentration of ${[S_0]}_0/{[R]}_0=0.2$ than for the higher ${[S_0]}_0/{[R]}_0=1$ and the lower ${[S_0]}_0/{[R]}_0=0.02$ value. The lower PS initial concentration ${[S_0]}_0/{[R]}_0=0.02$ is not high enough to produce the necessary amount of singlet oxygen ${}^1{\mathrm{O}}_2$ earlier, and the evolution of all the concentrations is slow. For the higher PS initial concentration ${[S_0]}_0/{[R]}_0=1$, a phenomenon known as self-shielding occurs. The PS concentration $[S_0]$ is so high that the light does not penetrate deep into the irradiated volume and the concentration of oxygen $[{}^3{\mathrm{O}}_2]$ is not high enough to support the process with photobleaching. Once the oxygen $[{}^3{\mathrm{O}}_2]$ is fully depleted, photobleaching stops and the successfully treated area $A_{\mathrm{success}}$ does not increase further. The photon density $\rho (\overrightarrow{r},120\mathrm{s})$, the concentration of PS $[S_0](\overrightarrow{r},120\mathrm{s})$, and the concentration of oxygen $[{}^3{\mathrm{O}}_2](\overrightarrow{r},120\mathrm{s})$ after 120 s are shown for this set of initial condions (${[S_0]}_0/{[R]}_0=1$) in Fig. 10. The oxygen concentration $[{}^3{\mathrm{O}}_2](\overrightarrow{r},120\mathrm{s})$ is reduced to 0 in the upper part of the dental pocket ${\mathrm{\Omega }}_{\mathrm{aPDT}}$, but the PS concentration $[S_0](\overrightarrow{r},120\mathrm{s})$ is not drastically reduced in this area. Only close to the border to air ${\mathrm{\Gamma }}_{\mathrm{air}}$ is the PS concentration $[S_0](\overrightarrow{r},120\mathrm{s})$ reduced to 0. This is enabled by the diffusion of oxygen from the air into the dental pocket ${\mathrm{\Omega }}_{\mathrm{aPDT}}$. In experiments, self-shielding of too high a PS concentration has been observed for cancer therapy.³³ The qualitative behavior of our model is supported by these experiments.

Fig. 10

Photon density $\rho (\overrightarrow{r},120\mathrm{s})$, the concentration of PS $[S_0](\overrightarrow{r},120\mathrm{s})$, and the concentration of oxygen ${[}^3{\mathrm{O}}_2](\overrightarrow{r},120\mathrm{s})$ after 120 s of irradiation. The PS initial condition is ${[S_0]}_0/{[R]}_0=1$. All the other computation parameters are given in Table 1.

The impact of coupling radiative transfer and the ongoing rate equations can be seen in Fig. 11. We compared the successfully treated area $A_{\mathrm{success}}/A_0$ of our model with continuously updated optical parameters to a model with constant optical parameters. The first Monte Carlo picture served as input for the spatially distributed photon density in the case of constant optical properties. The time when the models show a successfully treated area $A_{\mathrm{success}}/A_0>0$ different from zero is the same for both models. During this period, the concentrations change slightly and therefore an update of the optical properties in the irradiated area and thus of the photon density distribution is not necessary. Beyond this point in time, the successfully treated area $A_{\mathrm{success}}/A_0$ increases faster for our model with updated optical properties than for the calculation with the constant light distribution. For our model, the successfully treated area $A_{\mathrm{success}}/A_0$ increases nearly linearly in time for the first 38 s. Afterward, the laser radiation reaches the gingiva ${\mathrm{\Omega }}_{\mathrm{ging}}$ below the dental pocket ${\mathrm{\Omega }}_{\mathrm{aPDT}}$. A high percentage of the photons reaching this area ${\mathrm{\Omega }}_{\mathrm{ging}}$ is scattered back into the dental pocket ${\mathrm{\Omega }}_{\mathrm{aPDT}}$ and accelerates the aPDT there. Thereby, the steep increase after approximately 40 s in Fig. 11 is explained. After 47 s of irradiation, aPDT was successful in the whole dental pocket ${\mathrm{\Omega }}_{\mathrm{aPDT}}$.

Fig. 11

Ratio of successfully treated and total area $A_{\mathrm{success}}/A_0$ evolving in time for our model presented here with continuously updated optical properties and for the reaction diffusion model with a temporal constant photon density $\rho (\overrightarrow{r})$. The result of the first Monte Carlo run is used as this constant photon density $\rho (\overrightarrow{r})$. The computation parameters are given in Table 1.

In order to give an insight into the evolution of the other quantities, the photon density $\rho$, the survival fraction $B$, and the concentrations of the three PS states, the two oxygen states, and the receptors are shown in Fig. 12 after 40 s of irradiation from the top. (Movie 1 for 50 s of irradiation in online version).

Fig. 12

Simulation of aPDT in a dental pocket with irradiation from top for 40 s. The photon density $\rho$, the survival fraction $B$, and the concentrations of the three PS states, the two oxygen states, and the receptors are displayed. The parameters are taken from Table 1. (Movie 1 (MOV 2.26 MB) [URL: http://dx.doi.org/10.1117/1.JBO.19.7.071411.1] for 50 s of irradiation in online version).

3.3.

Treatment Protocols

In this paper, we present the simulated therapy progress for two paradigmatic treatment protocols. In Figs. 13 and 14, the therapy progress for irradiation from the top of the dental pocket ${\mathrm{\Gamma }}_{\mathrm{abs},1}$ can be seen. The used parameters are given in Table 1. In Fig. 13, the photon density $\rho (\overrightarrow{r},t)$ in the calculation area is shown at the beginning of irradiation (left) and after 40 s of irradiation (right). The survival fraction of bacteria $B(\overrightarrow{r},40\mathrm{s})$ after 40 s of irradiation is shown in Fig. 14.

Fig. 13

Simulation of aPDT in a dental pocket with irradiation from top. (a) photon density $\rho (\overrightarrow{r},0\mathrm{s})$ after 0 s of irradiation. (b) photon density $\rho (\overrightarrow{r},40\mathrm{s})$ after 40 s of irradiation. The parameters are taken from Table 1.

Fig. 14

Simulation of aPDT in a dental pocket with irradiation from top. Survival fraction $B(\overrightarrow{r},40\mathrm{s})$ after 40 s of irradiation. The parameters are taken from Table 1.

Simulations of trans gingival irradiation from the side via ${\mathrm{\Gamma }}_{\mathrm{abs},2}$ are performed with a wider beam radius of $w_0=4\mathrm{mm}$. The rest of the set of parameters is given in Table 1. The results for the photon density $\rho (\overrightarrow{r},t)$ after 0 s (left) and 80 s (right) are shown in Fig. 15. As we increased the laser radius compared to the radiation from the top via ${\mathrm{\Gamma }}_{\mathrm{abs},1}$, but left the maximum intensity the same, the absorbed laser power is increased for the irradiation from the side. Nevertheless, after 40 s and even after 80 s of irradiation, the dental pocket ${\mathrm{\Omega }}_{\mathrm{aPDT}}$ is not totally treated successfully. In Fig. 16, the survival fraction $B(\overrightarrow{r},80\mathrm{s})$ is shown after 80 s of irradiation. By this time, some area of the dental pocket has been treated successfully. Obviously, the protocol for the transgingival treatment has to be planned carefully under the criteria that gingiva is not hurt by the laser irradiation, but the dental pocket is treated completely.

Fig. 15

Simulation of aPDT in a dental pocket with trans gingival irradiation from left side. (a) photon density $\rho (\overrightarrow{r},0\mathrm{s})$ after 0 s of irradiation. (b) photon density $\rho (\overrightarrow{r},80\mathrm{s})$ after 80 s of irradiation. The parameters are taken from Table 1.

Fig. 16

Simulation of aPDT in a dental pocket with transgingival irradiation from left side. Survival fraction $B(\overrightarrow{r},80\mathrm{s})$ of bacteria after 80 s of irradiation. The laser radius is $w_0=4\mathrm{mm}$. The other parameters are taken from Table 1.