Modeling of Fundus Laser Exposure for Estimating Safe Laser Coagulation Parameters in the Treatment of Diabetic Retinopathy

A personalized medical approach can make diabetic retinopathy treatment more effective. To select effective methods of treatment, deep analysis and diagnostic data of a patient’s fundus are required. For this purpose, flat optical coherence tomography images are used to restore the three-dimensional structure of the fundus. Heat propagation through this structure is simulated via numerical methods. The article proposes algorithms for smooth segmentation of the retina for 3D model reconstruction and mathematical modeling of laser exposure while considering various parameters. The experiment was based on a two-fold improvement in the number of intervals and the calculation of the root mean square deviation between the modeled temperature values and the corresponding coordinates shown for the convergence of the integro-interpolation method (balance method). By doubling the number of intervals for a specific spatial or temporal coordinate, a decrease in the root mean square deviation takes place between the simulated temperature values by a factor of 1.7–5.9. This modeling allows us to estimate the basic parameters required for the actual practice of diabetic retinopathy treatment while optimizing for efficiency and safety. Mathematical modeling is used to estimate retina heating caused by the spread of heat from the vascular layer, where the temperature rose to 45 ◦C in 0.2 ms. It was identified that the formation of two coagulates is possible when they are located at least 180 μm from each other. Moreover, the distance can be reduced to 160 μm with a 15 ms delay between imaging.


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. 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-Al 2 O 3 /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.

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.
were investigated: an isothermal wall and isothermal flow. On the basis of numer methods, the authors considered systems of nonlinear ODEs which were solved using sequential relaxation method and Simpson's rule. They investigated different mesh mensions using numerical methods. The increase in the dimensions allowed the auth to gradually obtain more accurate solutions.
Despite the widespread use of numerical methods, including the integro-interp tion method, to solve various applied problems, their application to 3D modeling of f dus laser exposure involves significant computational difficulties. The development o gorithms that efficiently use modern High Performance Computing resources allowed to make significant progress in resolving this problem.
In the article, we discuss algorithms and methods for analyzing OCT images of fundus and propose algorithms for the mathematical modeling of laser exposure to fundus to estimate the appropriate parameters for safe and effective treatment. Moreo this article is the first to study the convergence of the integro-interpolation method for problem of mathematical modeling of laser exposure to the fundus. Thus, we develo a new 3D fundus structure model based on processing of OCT images and using app imating functions to describe the boundaries of the retina. One of the most important tors determining the safety of treatment is the distance between coagulates and the de between laser shots. In this case, one of the options for evaluating safe parameters is m ematical modeling. Indeed, simulating heat propagation at various values with the no parameters allowed us to estimate the basic parameters required for the actual practic diabetic retinopathy treatment while optimizing for efficiency and safety.

Research Material
A diabetic macular edema is one of the most unfavorable consequences of diab retinopathy, and can cause blindness. Figure 1 shows diagnostic images of a healthy f dus and a fundus damaged by macular edema. The fundus has a three-dimensional structure. For this reason, flat image anal may not always be effective. Currently, the OCT scanning of the fundus provides secti of the retinal image obtained in the oXZ plane. As a result, 85 cross-sectional images issued for different positions of the Y coordinate. Thus, the three-dimensional structur the fundus can be reconstructed using the data of a sequence of images, given the exa value in each OCT image. Figure 2 shows one of the sections obtained by OCT registra and also schematically shows the laser contacting the fundus. 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 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.
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.

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

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]. Mathematics 2021, 9, x FOR PEER REVIEW 6 of 16 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].
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 three-dimensional 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]: Step 1. Constructing a halftone image.
Step 2. Detection of the contour lines of the retina.
Step 3. Selection of a group of points located on contour lines.
Step 4. Approximation of contour lines through selected points.
Step 5. Construction of the image with selected layers based on smoothed contour lines.
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 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]: Step 1. Constructing a halftone image.
Step 2. Detection of the contour lines of the retina.
Step 3. Selection of a group of points located on contour lines.
Step 4. Approximation of contour lines through selected points.
Step 5. Construction of the image with selected layers based on smoothed contour lines.
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. 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. 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 two-dimensional OCT images makes it possible to reconstruct three-dimensional models of the fundus.

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): 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: Ir t x y z T C where c T is the temperature at the time of the laser shot, ( , , ) x y z β β = is the environment absorption function, ( , , ) v v C C x y z = is a function of the environment volumetric heat capacity at a fixed timestamp, x y z where the laser coagulation was initialized, and t Δ is the laser exposure duration. 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.

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): 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: where T c is the temperature at the time of the laser shot, β = β(x, y, z) is the environment absorption function, C v = C v (x, y, z) is a function of the environment volumetric heat capacity at a fixed timestamp, r = (x − x 0 ) 2 + (y − y 0 ) 2 + (z − z 0 ) 2 is the distance to the point (x 0 , y 0 , z 0 ) where the laser coagulation was initialized, and ∆t 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): where T = T(x, y, z, t) is temperature distribution depending on the spatial and time coordinates, k = k(x, y, z, T) is a function of thermal conductivity of the environment in space and time, C v = C v (x, y, z, T) 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, grad xyz is an operator for calculating the gradient of a function by coordinates x, y, z, and T 0 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 ψ(x, y, z) 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 C v 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 T = T + T 0 , where T 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 .
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 , y 0 , z 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 (T sr = 38 − 40 • 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).
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.
The iterative solutions for these tasks are given as follows: First, a solution to Equation (6) is sought on the time interval [t k , t k+1 ]. In this case, the initial condition of Equation (6) at moment t k 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 t k is nothing more than the result of solving Equation (6) at the time t k+1 or W t k+1 . It should be noted that for the splitting method the solution to the main problem, namely, temperature distribution at time t k+1 , can be represented as a solution to task (7), i.e., T| t=t k+1 ≈ V| t=t k+1 .
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++.

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. and (7). In this case, the resulting two-dimensional task can also be decomposed into onedimensional 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++.

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. 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 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][25][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 2N 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. 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 0 and z 0 , and research was carried out for different x 0 values. Figure 6 shows the results of such simulations for three different x 0 scenarios. It should be noted that the assessment was carried out with absolute parameters characterizing the displacement distance x 0 from the left edge, where x 0 = 0. 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 0 y and 0 z , and research was carried out for different 0 x values. Figure 6 shows the results of such simulations for three different 0 x scenarios. It should be noted that the assessment was carried out with absolute parameters characterizing the displacement distance 0 x from the left edge, (a) (b) (c) 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. 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 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.  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 7 shows that the retina heats up rather quickly. It should be noted 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 vascula the temperature on the retina could reach 44.9 °C, after which there would be a s crease in temperature with a limit trending to the normal temperature of the tiss clear that if one more shot is fired during these 0.24 ms, then the peak temperatu even exceed 44.9 °C. Based on this model, the required delay time between shots estimated.
The resulting model allowed us to simulate two shots in the area of the ep layer. Figure 8 shows an example of the implementation of the model. The coo along the oY axe was fixed, and the points of the shot were specified with differen along the oX and oZ axes. It should be noted that the different colors correspond to temperatures, so the figure can be interpreted in the same way in K and °C. Analysis of Figure 8 shows that as a result of each shot, heat spread from the lial layer to the retinal layer. A short time after exposure, the temperature of the re 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. 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.

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

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.