**1. Introduction**

The assessment of both the result of a service rendered to a person and the way it is provided is becoming extremely relevant; however, this approach seems far from being well developed, though modern computing systems have allowed the customization of unique patient treatment parameters by taking into account patient-specific features. This helps to increase treatment efficiency and patient opinion regarding the course of treatment. The application of personalized methods and means in treatment, both in diagnosis and monitoring, ensures the maximum treatment effect. In a broad sense, the concept of personalized medicine can be applied to specific diseases. In particular, taking into account the structure of a patient's fundus and the parameters of laser exposure will help to increase the number of successful operations in the future.

**Citation:** Shirokanev, A.; Ilyasova, N.; Andriyanov, N.; Zamytskiy, E.; Zolotarev, A.; Kirsh, D. Modeling of Fundus Laser Exposure for Estimating Safe Laser Coagulation Parameters in the Treatment of Diabetic Retinopathy. *Mathematics* **2021**, *9*, 967. https://doi.org/ 10.3390/math9090967

Academic Editor: Lucas Jódar

Received: 2 April 2021 Accepted: 24 April 2021 Published: 26 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Laser coagulation of the retina in the treatment of a diabetic macular edema is widespread today [1]. It involves multiple exposures to the edema to form coagulates [2]; however, coagulate deposition may have negative effects. When choosing unsuitable laser exposure parameters and locations for coagulates, the coagulation can lead to retinal damage and complete blindness [3]. This is because the use of a laser with increasing intensity leads to strong heating of a large area beyond the boundaries of the diabetic macular edema area, which, in turn, leads to the formation of a coagulum in the healthy retina during its critical heating. The maximum therapeutic effect becomes significantly complicated with an increase in the distance between coagulates, as well as under the conditions of lasers acting on blood vessels, retinal hemorrhages, and "solid" exudates [4]. These drawbacks in treatment, one way or another, often lead to irreversible decreases in vision. The best treatment can be obtained when coagulates are located at the same distance from each other and do not extend beyond the edema area, where the anatomical structures of the retina are not affected by the laser exposure.

Three modes provide the formation of coagulants under laser action on the fundus [5]: manual, semiautomatic, and navigation. Manual mode implies pointed laser shots which form respective coagulates in units. Semi-automatic mode denotes a series of shots leading to the formation of a group of coagulates in accordance with preset template areas. Thus, the locations of coagulates significantly affect treatment effectiveness; however, it is clear that a universal plan for their arrangement does not exist. The most effective one to date depends on the specific fundus structure and the location of the diabetic macular edema. The effectiveness of the placement is estimated by the locations of the centers of the coagulates formed. The most advanced means for arranging coagulates today is the Navilas unit, which provides high therapeutic efficiency. Operation in navigation mode requires the use of a preliminary coagulate arrangement plan based on a fundus image [6].

This approach does not solve at least two problems. The first is how to build a preliminary plan for coagulates, and the second is defining the optimal power and time of exposure for the ocular fundus. To solve the first problem, algorithms must be developed for the analysis of optical coherence tomography (OCT) images to isolate sensitive areas in the fundus, detect diabetic macular edema areas, and then choose the locations for coagulates [7]. For this purpose, various algorithms for pattern recognition, object detection, and image segmentation [8–11] can be used. To cope with the second problem, mathematical modeling has been proposed to specify the laser exposure parameters [12]. Indeed, on the one hand, it is necessary to achieve a predetermined temperature in the area of a diabetic macular edema for the formation of a coagulum. On the other hand, it is required to ensure a permissible temperature in the healthy area of the retina and other elements of the fundus to avoid negative effects. Consequently, the use of a mathematical model for heat propagation along the fundus should lead to an ideal temperature distribution in the fundus in a certain period, depending on the locations of coagulates, the power and duration of laser exposure, the duration of pauses between shots, etc.

Today, special attention in research is directed to the development of algorithms for processing fundus images [13–17]. Most importantly, the analysis of patient images using computer vision techniques in the treatment of diabetic retinopathy makes the process personalized and therefore more efficient and safer. We briefly consider the current state of research in the field of laser therapy.

In Reference [18], the estimation of parameters for laser coagulation was considered, but for the treatment of varicose veins. The authors noted unsafe treatments at high power radiation in the range of 8–20 W with a wavelength of 810–980 nm and 5–15 W for electromagnetic waves with a 1470–1550 nm wavelength. The study simulated the laser action of a solid-state laser. The advantages of the approach in comparison with real experiments was noted. Based on the obtained temperature values, the range of the permissible exposure power was identified. The lower and upper power levels were set accordingly. At the low level, the required effect was achieved, i.e., thrombus formation. The upper level had

power peaks, which caused irreversible damage to veins. The disadvantage of this work is that the model describes the operation of only one laser of a specific type and brand.

A detailed review of various lasers was presented in Reference [19]. The authors carried out an in-depth analysis of studies of various characteristics of lasers around the world. They mainly compared argon or diode lasers, as well as lasers based on panretinal photocoagulation (PRP) technology. The analysis showed that argon lasers, which have been considered the standard in the treatment of diabetic retinopathy for a long time, today no longer always provide the best effects.

The work in Reference [20] was conducted to assess anatomical and functional outcomes for patients undergoing treatment for diabetic retinopathy with a 532 nm laser PRP method (Supra Scan® Quantel Medical) and laser pinpoint coagulation at a wavelength of 532 nm (PASCAL® Topcon). The study was carried out with 48 patients and aimed to identify the correlation of their comfort with the laser exposure parameters, thereby showing the advantage of a multipoint laser.

A study of the thickness of the layer of retinal nerve fibers in the context of laser action was reported in Reference [21]. The article discusses lasers based on PRP technology. The effectiveness of red laser radiation treatment was compared to that of green laser treatment. The red laser showed an increase in the thickness of the retinal nerve fiber layer by 3–10 microns, which is 1.5 times higher than the green laser facilitated. The study also showed that the increase in thickness usually does not coincide with age or the number of burns.

In addition, Reference [22] showed that when comparing pain scores for patients with diabetic retinopathy, new systems, such as the novel navigated laser (NNL) system, have an advantage over conventional PRP-based systems.

The work of Reference [23] concerns laser action efficiency in terms of promoting the therapeutic treatment of diabetic retinopathy. Single-spot laser (SSL) and pre-stabilized laser (PSL) methods were compared, where the latter showed the best performance in terms of pain and effectiveness of treatment.

The analysis of the literature has shown that a limited number of works have been devoted to the problem of the mathematical modeling of fundus laser exposure. Nevertheless, this problem can now be solved via the use of numerical methods. In Reference [24], numerical methods based on an implicit difference scheme were considered. Such methods are often used to obtain fuzzy wave equations. The main idea of the work was the formulation of a new space of coordinates with respect to time using implicit schemes and the theory of fuzzy sets. The use of fuzzy sets allowed the authors to perform fuzzy analysis of the resulting equations and ensured function convergence. The use of numerical methods allowed the authors to reduce the time required for constructing wave equations with minimal losses in accuracy when describing real processes. The confirmation of the convergence of the methods made it possible to assess the possibility of using such mathematical models with a given error.

Indirectly, thermal energy distribution and transfer mechanics were studied in Reference [25], where a model of the flow of dusty nanofluid Cu-Al2O3/water was developed. Using numerical methods, the authors obtained ordinary differential equations (ODEs) that described the physical process of such a flow. In this case, double solutions were obtained, which were investigated for convergence and after which only one solution was chosen as a reliable one. On that basis, various parameters of the process were evaluated. For example, it has been found that nanosized particles (Cu) have a significant effect on heat transfer.

Finally, Reference [26] also studied heat propagation; however, the authors investigated the dynamics of micropolar fluids using numerical methods. Two edge conditions were investigated: an isothermal wall and isothermal flow. On the basis of numerical methods, the authors considered systems of nonlinear ODEs which were solved using the sequential relaxation method and Simpson's rule. They investigated different mesh

dimensions using numerical methods. The increase in the dimensions allowed the authors to gradually obtain more accurate solutions.

Despite the widespread use of numerical methods, including the integro-interpolation method, to solve various applied problems, their application to 3D modeling of fundus laser exposure involves significant computational difficulties. The development of algorithms that efficiently use modern High Performance Computing resources allowed us to make significant progress in resolving this problem.

In the article, we discuss algorithms and methods for analyzing OCT images of the fundus and propose algorithms for the mathematical modeling of laser exposure to the fundus to estimate the appropriate parameters for safe and effective treatment. Moreover, this article is the first to study the convergence of the integro-interpolation method for the problem of mathematical modeling of laser exposure to the fundus. Thus, we developed a new 3D fundus structure model based on processing of OCT images and using approximating functions to describe the boundaries of the retina. One of the most important factors determining the safety of treatment is the distance between coagulates and the delay between laser shots. In this case, one of the options for evaluating safe parameters is mathematical modeling. Indeed, simulating heat propagation at various values with the noted parameters allowed us to estimate the basic parameters required for the actual practice of diabetic retinopathy treatment while optimizing for efficiency and safety.

### **2. Materials and Methods**

### *2.1. Research Material*

A diabetic macular edema is one of the most unfavorable consequences of diabetic retinopathy, and can cause blindness. Figure 1 shows diagnostic images of a healthy fundus and a fundus damaged by macular edema.

**Figure 1.** Fundus images: (**a**) healthy fundus; (**b**) fundus damaged by diabetic macular edema.

The fundus has a three-dimensional structure. For this reason, flat image analysis may not always be effective. Currently, the OCT scanning of the fundus provides sections of the retinal image obtained in the oXZ plane. As a result, 85 cross-sectional images are issued for different positions of the Y coordinate. Thus, the three-dimensional structure of the fundus can be reconstructed using the data of a sequence of images, given the exact Y value in each OCT image. Figure 2 shows one of the sections obtained by OCT registration and also schematically shows the laser contacting the fundus.

The edge condition is an area where the laser action has no effect on the vascular layer. The presence and positioning of such an area is extremely important for laser exposure modeling. As the analysis of the subject area has shown, the identification of new methods of effective treatment for diabetic retinopathy is possible when a sufficiently large number of fundus images can be analyzed. Additional features can be extracted from the analysis

of OCT images. This requires a lengthy and thorough analysis of a sample of images with the involvement of medical experts.

**Figure 2.** The structure of the retina as it pertains to laser exposure: 1—edge conditions, 2—vitreous, 3—retina, 4—vascular layer, 5—laser radiation.

The collection of a large database for the assessment of laser exposure effects on the fundus seems to be a great challenge. At present, the only way to assess the effective parameters of laser exposure is to carry out laser coagulation and fix the values under various conditions. This approach requires an excessively large number of patients as samples, since the physician cannot test different laser exposure parameters on real patients to study the effectiveness of the laser coagulation. A small sample of treated patients may not contain enough information to identify safe parameters. One approach is recommended that would require a small patient base containing heterogeneous fundus structures. As a result, the optimal material for research is a mathematical model that describes the spread of heat along the fundus in three-dimensional space, depending on the laser power, the duration of the exposure, etc. Having defined the temperature inside the fundus under normal conditions and the temperature of the formation of coagulates, the optimal parameters of laser exposure can be determined, including the target area. This makes it possible to explore various changes in the action of an object by replacing it with a mathematical model.

#### *2.2. OCT Image Analysis and Reconstruction of 3D Fundus Structure Model*

To assess the parameters of laser exposure, the preliminary analysis of OCT images may be required. Greater accuracy in determining the main layers on the cross-sections of the fundus images will result in a more suitable reconstructed model of the fundus surface. The reconstructed 3D model is needed to estimate the heat propagation.

At present, there are no universal methods for registering a fundus that allow obtaining undistorted images of parts of the eye. For various reasons, the recorded areas of the eye's surface in the original image may contain additional interference (Figure 3a).

Consequently, the first task is to highlight different layers in the noisy fundus image. The preprocessing stage may include the filtering procedures. For example, as mentioned earlier, median filtering ensures the elimination of strong impulse noise by brightness replacement with the median value in a certain neighborhood. If the noise is additively distributed over all pixels in the form of additive Gaussian white noise (AGWN), then it is possible to use recurrent Kalman filtering [27].

**Figure 3.** The result of the optical coherence tomography (OCT) fundus image processing: (**a**) original OCT image against the background of noise; (**b**) preliminary segmentation of the retinal layer; (**c**) fundus image model based on parametric functions.

Thus, preprocessing does not solve a specific applied problem but improves the quality of the solution obtained in the course of the main image processing. Based on the image shown in Figure 3a, the segmentation of three layers is required: the vitreous, retinal, and epithelial (vascular) layers; however, using the a priori knowledge that the retina always separates the other two layers, the task is simplified to binary segmentation, where the retinal layer will be highlighted. With segmentation, it is important to take into account that the retinal layer must be evenly filled and distributed over the entire width of the image. In the image, every pixel inside the retina should have zero brightness, and every pixel outside the retina should have maximum brightness. At the first stage, the manual selection of borders and filling of the corresponding area can be performed. Figure 3b shows the result of such processing as a binary image.

The analysis of Figure 3b shows a rather complex structure of the boundaries of the retinal layer, which can lead to significant difficulties in the reconstruction of a threedimensional model. Therefore, the next processing step is the generation of a binary image, for which the boundary lines of the retinal layer can be approximated using mathematical functions and have a smooth appearance. Continuous parametric functions are used to describe the upper and lower boundaries of the retinal layer. Figure 3c shows the resulting image after approximation.

As can be seen from Figure 3c, in addition to smoothing, three areas of the fundus are segmented separately in the final modeled image.

Thus, after preliminary processing, the considered algorithm will consist of the following five steps [28]:


The sequential processing of all OCT images for the fundus allows one to obtain a set of approximating functions and use them to restore the three-dimensional structure of the fundus from OCT images. Figure 4 shows an example of such reconstruction. The processing of OCT images and reconstruction of a 3D model was performed using an implementation of the proposed algorithm in MATLAB.

**Figure 4.** Fundus model reconstructed from OCT images.

Using 3D fundus models, based on the modeling of laser exposure, it is possible to analyze to what temperature this or that point of the fundus will be heated during laser exposure with certain parameters. This allows the prediction of the formation of coagulates in accordance with the chosen plan and to assess the likelihood of unwanted side effects, including heating the retinal layer to a critical temperature.

The analysis of Figure 4 shows that the use of approximating functions for twodimensional OCT images makes it possible to reconstruct three-dimensional models of the fundus.

#### *2.3. Heat Propagation Modeling*

In the mathematical modeling of laser exposure, it is necessary to take into account the fact that energy is converted in such a model, i.e., a laser light pulse (light energy) leads to heating of the fundus (thermal energy). Indeed, the electromagnetic energy of the laser action is converted into heat when interacting with the vascular layer [29,30].

The intensity of light energy is described by Equation (1):

$$I(r) = \frac{P}{\pi a^2} e^{-\left(\frac{r}{a}\right)^2},\tag{1}$$

where *r* is the distance from the light source, *a* is the spot radius, and *P* is the light source's power.

Using Equation (2), the temperature distribution in three-dimensional space can be determined after the end of the laser exposure:

$$\Psi(x, y, z) = \frac{e^{-\frac{z}{\hbar}\beta(x, y, z^{\sharp})d\xi}}{\mathbb{C}\_{\upsilon}}\beta I(r)\Delta t \tag{2}$$

where *Tc* is the temperature at the time of the laser shot, *<sup>β</sup>* <sup>=</sup> *<sup>β</sup>*(*x*, *<sup>y</sup>*, *<sup>z</sup>*) is the environment absorption function, *Cv* <sup>=</sup> *Cv*(*x*, *<sup>y</sup>*, *<sup>z</sup>*) is a function of the environment volumetric heat capacity at a fixed timestamp, *<sup>r</sup>* = / (*<sup>x</sup>* <sup>−</sup> *<sup>x</sup>*0) 2 + (*<sup>y</sup>* <sup>−</sup> *<sup>y</sup>*0) 2 + (*<sup>z</sup>* <sup>−</sup> *<sup>z</sup>*0) <sup>2</sup> is the distance to the point (*x*0, *<sup>y</sup>*0, *<sup>z</sup>*0) where the laser coagulation was initialized, and <sup>Δ</sup>*<sup>t</sup>* is the laser exposure duration.

Considering very small values of Δ*t*, it is possible to neglect the diffraction of light [27]. Since this is indeed the case in practice, the model of the heat propagation after laser exposure can be rewritten in the form of Equation (3):

$$\begin{cases} \ c\_v \frac{\partial T}{\partial T} = \operatorname{div} \left[ k \cdot \operatorname{grad}\_{xyz}(T) \right], \\ \ T|\_{t=0} = \psi(x, y, z), \\ \ T|\_E = T\_{0\_\prime} \end{cases} \tag{3}$$

where *<sup>T</sup>* <sup>=</sup> *<sup>T</sup>*(*x*, *<sup>y</sup>*, *<sup>z</sup>*, *<sup>t</sup>*) is temperature distribution depending on the spatial and time coordinates, *<sup>k</sup>* <sup>=</sup> *<sup>k</sup>*(*x*, *<sup>y</sup>*, *<sup>z</sup>*, *<sup>T</sup>*) is a function of thermal conductivity of the environment in space and time, *Cv* <sup>=</sup> *Cv*(*x*, *<sup>y</sup>*, *<sup>z</sup>*, *<sup>T</sup>*) is a function of the environment volumetric heat capacity which changes during heating or cooling, *E* is the edge laser exposure influence on temperature, *div* is a vector field divergence operator, *gradxyz* is an operator for calculating the gradient of a function by coordinates *x*, *y*, *z*, and *T*<sup>0</sup> is the temperature at the edge region.

When using Equation (3), it is necessary that the region of determination of the thermal field in space is large enough to ensure that, under laser action, the propagation of heat only occurs up to the edge through which the laser passes. Two methods can be used to meet this condition.

First, the function *<sup>ψ</sup>*(*x*, *<sup>y</sup>*, *<sup>z</sup>*) must be specified at the bearing edge, i.e., one that changes over time. In addition, it is necessary to perform a special linear replacement in order to reduce the task to fixed (or zero) edge conditions. Thus, the thermal conductivity will be described by an inhomogeneous differential equation. Solution retrieval is possible via expansion into a solution for the corresponding homogeneous equation with the initial edge and initial conditions. Next, a solution to the inhomogeneous equation is found for which a zero edge and initial conditions are provided. Consequently, the final solution will be found from an equation in the form of Equation (3).

Nevertheless, several difficulties usually arise when searching for such a replacement and solving an inhomogeneous equation. Therefore, in this work, it was decided to use the second method. This method is based on the fact that it is necessary to artificially expand the domain of definition, and then, based on the resulting extension, divide it into informative and non-informative areas. In this case, for the first region, the condition of a relatively negligible pulse duration is satisfied, which makes it possible to describe the temperature distribution at the initial moment of time in accordance with Equation (2). Determining the spread of heat in a non-informative area will not be difficult, since for solving this problem it will be sufficient to use a symmetric display; however, the already fixed temperature value will correspond to the border in the uninformative area. At the same time, when analyzing a mathematical model, it will be possible to use only data in the informative area which further simplifies the task.

Unfortunately, due to the dependence of the volumetric heat capacity and temperatureconductivity on temperature, the task will be nonlinear. Nevertheless, it is possible to assess the change in the shape of the retina based on the temperature values of its layers, i.e., by assessing the given layer and temperature. Equation (4) is similar to Equation (3) in which the functions of the volumetric heat capacity *Cv* and thermal diffusivity *k* depend only on spatial coordinates.

The analysis of Equation (4) shows that it is not difficult to obtain zero edge conditions. Indeed, by replacing *<sup>T</sup>* = *<sup>T</sup>*<sup>&</sup>gt; + *<sup>T</sup>*0, where *<sup>T</sup>*<sup>&</sup>gt; describes the simulated temperature in the form of the direct effect of laser exposure (i.e., how much a given point of the fundus has heated as a result of laser exposure). Absolute temperatures on the fundus surface *T* are defined as the sum of the heating temperature and the initial temperature *T*0.

$$\begin{cases} \mathbb{C}\_{v}(\mathbf{x}, \mathbf{y}, z) \frac{\partial T}{\partial t} = \operatorname{div} \left[ k(\mathbf{x}, \mathbf{y}, z) \cdot \operatorname{grad}\_{\mathbf{x}\mathbf{y}z} (T) \right], \\ \qquad T \Big|\_{t=0} = \frac{\epsilon^{-\beta(\mathbf{x}, \mathbf{y}, z)} \beta(\mathbf{x}, \mathbf{y}, z) I(r) \Delta t}{\mathbb{C}\_{v}(\mathbf{x}, \mathbf{y}, z)} + T\_{\epsilon}, \\ \qquad \qquad T|\_{E} = T\_{0}. \end{cases} \tag{4}$$

A normal tissue temperature (~36.5 ◦C) can be used as the starting temperature. The solution of Equation (4) makes it possible to obtain a model of the temperature change after the application of a laser with a given power at the point (*x*0, *<sup>y</sup>*0, *<sup>z</sup>*0) throughout the structure of the fundus. Then, by combining the resulting model with a three-dimensional fundus model, it is possible to predict the absolute temperature at each point of the reconstructed model. The temperature on the retina is estimated based on this alignment. If, at some point in the retinal layer, the value exceeds the critical temperature (*Tsr* <sup>=</sup> <sup>38</sup> <sup>−</sup> <sup>40</sup> ◦C), then treatment with the given parameters may be unsafe. Otherwise, the temperature in the area of the diabetic macular edema is estimated. If this temperature exceeds 39 ◦C, then a coagulum will appear and the treatment can be considered effective.

The solution of the task in an analytical form is not possible; however, the temperature distribution can be estimated using numerical methods [31]. For example, based on the splitting scheme, it is possible to go from a multidimensional task to a one-dimensional one. This allows the use of an implicit scheme without solving the linear relationship system. The application of a sweep method significantly speeds up the solution retrieval for this task.

To use the splitting method, it is necessary to rewrite Equation (4) in the form of Equation (5).

$$\begin{cases} \begin{array}{c} s\_{\mathbb{P}}(\mathbf{x},\boldsymbol{y},\boldsymbol{z}) \frac{\partial T}{\partial \boldsymbol{x}} = \frac{\partial}{\partial \boldsymbol{x}} (k(\boldsymbol{x},\boldsymbol{y},\boldsymbol{z}) \frac{\partial T}{\partial \boldsymbol{x}}) + \\ + \frac{\partial}{\partial \boldsymbol{y}} (k(\boldsymbol{x},\boldsymbol{y},\boldsymbol{z}) \frac{\partial T}{\partial \boldsymbol{y}}) + \frac{\partial}{\partial \boldsymbol{z}} (k(\boldsymbol{x},\boldsymbol{y},\boldsymbol{z}) \frac{\partial T}{\partial \boldsymbol{z}}), \\ \qquad T|\_{t=0} = \boldsymbol{\psi}(\boldsymbol{x},\boldsymbol{y},\boldsymbol{z}), \\ \qquad \qquad T|\_{E} = T\_{0}. \end{array} \tag{5}$$

The numerical solution of Equation (5) using the splitting method is described in a more convenient form below. The idea behind the method is that sampling with the same step must first be performed for a given time segment. After that, the solution is reduced to solving iterative Equations (6) and (7). At the first stage, the coordinate is split off along the oY axis, since this is the main axis of the OCT images.

$$\begin{cases} \ c\_v(\mathbf{x}, y, z) \frac{\partial \mathcal{W}}{\partial t} = \frac{\partial}{\partial y}(k(\mathbf{x}, y, z) \frac{\partial \mathcal{W}}{\partial y}),\\ \mathcal{W}|\_{t=t\_k} = T|\_{t=t\_{k'}} \\ \mathcal{W}|\_E = 0 \end{cases} \tag{6}$$

$$\begin{cases} s\_v(x, y, z) \frac{\partial V}{\partial t} = \frac{\partial}{\partial x}(k(x, y, z) \frac{\partial V}{\partial x}) + \\ \qquad + \frac{\partial}{\partial z}(k(x, y, z) \frac{\partial V}{\partial z})\_\prime \\ \qquad V|\_{t = t\_k} = W|\_{t = t\_{k'}} \\ \qquad V|\_E = 0. \end{cases} \tag{7}$$

The iterative solutions for these tasks are given as follows: First, a solution to Equation (6) is sought on the time interval [*tk*, *tk*+1]. In this case, the initial condition of Equation (6) at moment *tk* coincides with the solution to Equation (5) at the same time. Further, for Equation (7), a solution is sought under the condition that the initial condition at the moment *tk* is nothing more than the result of solving Equation (6) at the time *tk*+<sup>1</sup> or *Wtk*+<sup>1</sup> . It should be noted that for the splitting method the solution to the main problem, namely, temperature distribution at time *tk*+1, can be represented as a solution to task (7), i.e., *<sup>T</sup>*|*t*=*tk*+<sup>1</sup> <sup>≈</sup> *<sup>V</sup>*|*t*=*tk*+<sup>1</sup> .

It is necessary to understand that Equation (6) describes a set of tasks taking into account the dependence of the sought functions on all spatial coordinates. Similar dependencies are observed for Equation (7). The solution to Equation (6) can be obtained using the finite difference method based on an implicit difference scheme.

To summarize, the three-dimensional problem is reduced to solving Equations (6) and (7). In this case, the resulting two-dimensional task can also be decomposed into one-dimensional tasks or solved using numerical methods.

To carry out computational experiments of laser exposure to the fundus, the proposed algorithm for heat propagation modeling was implemented in C++.

#### **3. Results**

Several important results were obtained as a result of solving the laser exposure mathematical modeling task on the basis of numerical methods. First, the modeling of the laser action of the fundus was carried out while taking into account its three-dimensional structure. Figure 5 shows the result of this simulation. In addition, the simulation results can be used to estimate the parameters of laser exposure. By superimposing the obtained three-dimensional model of heat distribution at a given moment in time on the fundus model under the condition of normal tissue temperature, it is possible to estimate the absolute values of temperature at a given moment in time. Moreover, it is possible to estimate how quickly the temperature is restored to normal after the termination of the laser action. In Figure 5, some heat dissipation can be observed in the vicinity of the center point of the shot.

**Figure 5.** Three-dimensional model of laser exposure.

The use of numerical methods, generally speaking, can lead to different simulation results for different mesh sizes. Moreover, the larger the mesh size is, the longer the simulation will take. This gives rise to the problem of investigating the convergence of solutions.

The presence of discontinuities in the coefficients of thermal conductivity and volumetric heat capacity leads to the absence of the existence of derivatives at points of discontinuity. The integro-interpolation method (IIM) [31,32] aims at eliminating this problem, and consists of calculating integrals over the limits of the neighborhoods of the grid nodes for the main equation of the problem and the initial and boundary conditions. As a result of calculating the integrals, a difference scheme is formed which is considered to be resistant to discontinuities. It is very difficult to analytically show the effectiveness of the IIM in comparison with the finite difference method, as existing works have only been aimed at analyzing particular problems related to heat conduction [24–26]; however, there is a way to experimentally evaluate the convergence of the method. This method involves the numerical simulation of the coagulation process for two grids in which the first grid consist of *N* intervals and the second one consist of 2*N* intervals.

Since the modeling is performed with three spatial coordinates and one time coordinate, a convergence study was carried out on the basis of modeling at different sampling steps along one of the axes. At the same time, for all nodes available in two neighboring models, the temperature discrepancies were estimated and the standard deviation or root mean square (RMS) values were calculated. The results of the convergence study are shown in Tables 1–4. Here *I*, *J*, *K* are the numbers of intervals in spatial coordinates *x*, *y*, *z*, respectively, and *S* is number of intervals in time coordinate *t*.


**Table 1.** Convergence for different mesh sizes along the oX axis.

**Table 2.** Convergence for different mesh sizes along the oY axis.


**Table 3.** Convergence for different mesh sizes along the oZ axis.


**Table 4.** Convergence for different mesh sizes along the time axis.


Analysis of the results presented in Tables 1–4 shows that with an increase in the number of grid dimensions, the results of mathematical modeling converge; however, different convergence rates are set for different coordinates. For example, for the oX axis, an increase in the number of intervals from 60 to 480 (8 times) leads to a decrease in the standard deviation by 87 times—but 29 times for the oY axis and less than 36 times for the oZ axis. A halving in the number of steps along the spatial coordinates leads on average to decreases in the standard deviation by 2.9–5.7 times along the oX axis, 2.2–4.5 times along the oY axis, and 1.7–5.9 times along the oZ axis. A stable decrease in the standard deviation indicates the convergence of the IIM. Most discontinuities exist in the oZ direction and the convergence is slow there. In the presence of discontinuities, the standard deviation usually decreases by no more than 50%. There are discontinuities in the oX and oY directions which are not as pronounced as in the Z direction. In these directions, the standard deviation can decrease by more than 50%, which is shown in the tables above. When studying the influence of the dimensions when passing from 200 intervals to 3200 intervals, a decrease in the deviation by a factor of 11 was observed; however, in terms of absolute values, the results showed far lower values for spatial coordinates with a value of approximately 5 ◦C. Reducing the sampling time step by a factor of two led to a decrease in the standard deviation of no more than 50%. On average, as can be seen in Table 4, the standard deviation decreased by 1.8 times. Thus, the algorithm has good convergence; that is, an increase in the grid dimension size will lead to an improvement in the accuracy of the solution.

Heat propagation with various laser displacements was also simulated. In particular, the coordinates along the oY and oZ axes were fixed as *y*<sup>0</sup> and *z*0, and research was carried out for different *x*<sup>0</sup> values. Figure 6 shows the results of such simulations for three different *x*<sup>0</sup> scenarios. It should be noted that the assessment was carried out with absolute parameters characterizing the displacement distance *<sup>x</sup>*<sup>0</sup> from the left edge, where *<sup>x</sup>*<sup>0</sup> <sup>=</sup> 0.

**Figure 6.** Simulation of thermal propagation at various coordinates for the point of influence along the oX axis: (**a**) 200 μm; (**b**) 350 μm; (**c**) 500 μm.

Analysis of the simulation results presented in Figure 6 shows that the structure (shape) of the retina had little effect on the spread of heat. For the considered example, the thickness of the retinal layer was approximately the same for all displacements; however, there were slight differences in the spread of heat along the fundus under the influence of different displacements. This difference should be taken into account when forming several coagulates side by side in order to prevent reaching a critical temperature on the retinal layer. Figure 7 shows the results of modeling the temperature distribution on the retinal layer after laser exposure at different points in time.

**Figure 7.** Simulation of temperature versus time.

Figure 7 shows that the retina heats up rather quickly. It should be noted that the value at time zero is associated with laser action on the adjacent layer.

It can be seen that within 0.24 ms, due to the spread of heat from the vascular layer, the temperature on the retina could reach 44.9 ◦C, after which there would be a slow decrease in temperature with a limit trending to the normal temperature of the tissue. It is clear that if one more shot is fired during these 0.24 ms, then the peak temperature may even exceed 44.9 ◦C. Based on this model, the required delay time between shots can be estimated.

The resulting model allowed us to simulate two shots in the area of the epithelial layer. Figure 8 shows an example of the implementation of the model. The coordinate along the oY axe was fixed, and the points of the shot were specified with different offsets along the oX and oZ axes. It should be noted that the different colors correspond to heating temperatures, so the figure can be interpreted in the same way in K and ◦C.

**Figure 8.** Modeling of laser action for the formation of two coagulates: 1—vitreous layer, 2—retinal layer, 3—epithelial layer, 4—vascular layer.

Analysis of Figure 8 shows that as a result of each shot, heat spread from the epithelial layer to the retinal layer. A short time after exposure, the temperature of the retina can increase by 5 ◦C due to the spread of heat from one shot. The second shot, taking into account the proximity of the distance between the coagulates, will also contribute to the spread of heat in the local neighborhood. During treatment, it is impossible to allow a critical temperature on the retinal layer to be reached.

In this regard, in order to identify safe treatment parameters, heat propagation was simulated for various distances between the centers of coagulates and delay times between shots. In particular, the situations of laser exposure with two shots with different displacements and delay times were simulated. Figure 9 shows the results of such a simulation at different coordinates along the oX axis. It should be noted that the temperature is calculated as the maximum temperature over the entire area as a result of simulating two laser shots.

Analysis of Figure 9 shows that a safe distance between coagulates is about 180 μm. Another option to ensure safe treatment is to increase the delay time between shots to 15 ms; the distance can be reduced by increasing the delay between shots. A distance of around 160 μm may be used with a delay greater than 15 ms.

Thus, the main results here are the mathematical models of laser exposure which make it possible to estimate the parameters of safe treatment. Notably, these models feature the property of convergence.

**Figure 9.** Dependence of the maximum temperature on the epithelial layer during the implementation of two shots on the delay between shots and the coordinates of the shots.

#### **4. Discussion**

The results obtained within the framework of the study can be taken into account for the treatment of diabetic retinopathy, and the models developed can be used in treatment to select parameters of laser exposure that will provide the maximum therapeutic effect. The use of the proposed fundus models and laser exposure parameters is another step towards personalized medical care. Indeed, the OCT images make it possible to restore the three-dimensional structure of the fundus where the laser exposure is simulated. The modeling of heating the retinal, vascular, and epithelial layers allows conclusions to be formed regarding safe shifts between laser shots and pause durations.

Many quantitative results have been obtained. In particular, at a given laser power, a safe pause of 15 ms was revealed with a distance between coagulates of 160 μm. In addition, a minimum safe distance between shots was 180 μm, which did not depend on the delay between shots. Convergence of the numerical methods used in modeling has been established with a decrease in the sampling step size by a factor of two, where the standard deviation stably decreases. A decrease in the step size by 50% along the spatial coordinates leads, on average, to decreases in the standard deviation of 2.9–5.7 times along the oX axis, 2.2–4.5 times along the oY axis, and 1.7–5.9 times along the oZ axis. A stable decrease in the standard deviation indicates the convergence of the integro-interpolation method. This becomes especially important when arranging coagulates not as a single formation, but as a group. Moreover, on the basis of mathematical modeling, it is potentially possible to adjust the plans of coagulates in order to maximize the therapeutic effects. In the future, we plan to study the placement of a group of coagulates, as well as the performance of the proposed algorithms.

**Author Contributions:** Methodology, A.S., E.Z. and A.Z.; writing—original draft, N.A., A.S. and N.I.; writing—review and editing, N.A. and D.K.; supervision, N.I.; funding acquisition, N.I. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by RFBR (grant number 19-29-01135); Ministry of Science and Higher Education of the Russian Federation within (Government project).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.
