Fractal characterization of retinal microvascular network morphology during diabetic retinopathy progression

The study aimed to characterize morphological changes of the retinal microvascular network during the progression of diabetic retinopathy.


| INTRODUC TI ON
According to the World Health Organization global report on DM from 2016, an estimated number of adults living with DM increased from 108 million in 1980 to 422 million 2014 and the global agestandardized prevalence of DM increased from 4.7% to 8.5% in the same period. 1 In accordance with this, the prevalence of DM has been increasing over time in the United States, especially in the adults older than 65. 2 In most of the world, life expectancy has doubled in the past century, 3 so the burden of diabetes is expected to rise significantly as the number of individuals with DM increases.
Since the treatment of DM complications can cause 11-fold increase in direct medical cost of type 2 diabetes, 4 early diagnosis and timely treatment of DM complications are extremely important. DR is a microvascular complication of diabetes and one of the leading causes of loss of vision in middle income and industrialized countries in the world. 5 Therefore, the American Diabetic Association currently recommends referral of individuals with type 1 diabetes (DM1) to their first ophthalmologic exam 5 years after the initial diagnosis of diabetes and with type 2 diabetes (DM2) at the time of initial diagnosis. 6 This recommendation is based on the fact that at the time of the first diagnosis of diabetes, more than 30% of patients already have microvascular complications, but more than 40% of them have no associated symptoms. 7,8 Fundoscopic examination and visualization of retinal vascular network can give us not only the information related to the presence and severity of diabetic retinopathy, but also on the condition of the microvascular network in the entire organism, and potential presence of other microvascular complications of DM. 9 Numerous studies have confirmed that microvascular disease represents an independent risk factor for development of cardiovascular complications in DM1 and DM2. [9][10][11] Moreover, Exalto et al have shown that individuals with severe diabetic retinopathy have significantly increased risk for the development of dementia. 12 Taken together, these reports show that the assessment of the retinal microvasculature in patients with diabetes can be important not only for the prevention and treatment of blindness, but also for the timely prevention and treatment of other microvascular complications of DM, as well as cardiovascular complications and possibly dementia. Therefore, there is a need for the development of a simple, noninvasive method for DR screening and DR severity staging in a large population in the primary care setting. Inclusion of family physicians in this process of DR screening focused on early detection of DR, and prevention of blindness has already been shown to be very effective. 13 One of the methods that could be used as a screening tool is analysis of 2-dimensional color images of retina captured by digital fundus cameras. These cameras are affordable, becoming smaller in size and portable, and are able to capture high resolution retinal images of exceptional quality. 14,15 The complexity of the microvascular network is affected by DR. The studies using fractal analysis to explore this phenomenon report conflicting results: while some researchers report increased, the others report decreased microvascular complexity in F I G U R E 1 Image processing work flow. Two-dimensional color retinal images of patients with diabetes mellitus from the publicly available DIARETDB1 and STARE databases were used in the study. 14,15,25,26 STEP 1: The images were segmented by the automatic method of blood vessel segmentation that originated from the RETINAL 2020 project. 28 STEP 2: The fine and the coarse filtering settings and thresholding were used on each image to generate one binarized image with the low, and the second image with the high number of vessel branching generations. STEP 3: vascular network morphology in different stages of diabetic retinopathy progression was compared by fractal analysis parametrization such as box/count dimension, lacunarity, and multifractals in the set of images with low number of branching generations, as well as in the set with a high number of branching generations DR. [16][17][18] Morphological features typically found in DR such as arteriolar pruning and neovascularization 19 can affect complexity of microvascular network morphology in the opposing ways. While arteriolar pruning could decrease complexity, neovascularization that is present only in the advanced stages of DR could increase complexity of the microvascular network. Therefore, the observed inconsistencies in these reports may be related to different stages of DR being studied. Different resolution of details in the images of vascular network used in these studies could be the cause of these disagreements as well. We show that the results of fractal analysis are indeed affected by the number of vessel branching generations included in the analysis. More importantly, the study describes a simple approach that could be used to detect severe stages of DR even in primary care settings. This could improve the stratification of patients with DR according to the risk for development of other diabetic microvascular complications and allow for their timely referral to appropriate specialists, improvement of treatment outcomes, as well as the reduction of medical costs.  Since the severity of DR in each image was not specified, three of our medical experts independently determined the stage of DR for each image in the DIARETDB1 database by following the Early Treatment Diabetic Retinopathy Study Research Group recommendations. 21 The group of medical experts consisted of one ophthalmologist and two medical doctors-general practitioners. In order to make sure they used consistent grading criteria, all experts were required to complete the online self-directed diabetic retinopathy grading course and successfully pass the competency-based exam. 22 The retinal images were initially staged into five groups: normal looking with no signs of DR (stage 0), mild (stage 1), moderate (stage 2), severe non-proliferative (stage 3), and proliferative DR (stage 4).

| Retinal images and diabetic retinopathy staging
In cases where there were slight differences in staging among the medical experts, the stage was determined based on the majority opinion. Since our study focused on a subset of only 20 samples that were adequately illuminated, the images were finally grouped into three groups according to severity of DR (Table 1): normal looking (five images with no signs of DR-stage 0), moderate DR (six images with signs of moderate DR-stage 2), and severe DR (five images with severe non-proliferative DR-stage 3 and 4 images with signs of severe proliferative DR-stage 4). In this group of 20 images, there were no images corresponding mild DR-stage 1. The overall agreement among the graders was 80%, with a good value for free-marginal kappa of 0.70 (95%CI = 0.49-0.91). This inter-rater variability in our study was comparable to other studies using the same method of measurement of inter-rater variability. 23,24 Among the three groups of images, there was no significant difference in the estimated light exposure and in the estimated tonal range (data not shown).
In addition, the raw retinal images, from the publicly available Structured Analysis of Retina (STARE) database, were used in the study. 25,26 The STARE database contains 402 annotated color images captured at 35-degree field of view with 700 × 605 pixel resolution. Out of these, we used a subset of images that were either normal looking or associated with the diagnosis of proliferative diabetic retinopathy (PDR) and that were of adequate quality as previously described. 27 In this subset, ten images did not have any pathological changes, while 9 were associated with the diagnosis of PDR. Some of the images were centered on the macula, but some were centered on the optic disk. All 19 images belonging to the STARE database were adequately illuminated, and the mean estimated exposure values ranged from 101.7 to 176.9. The cumulative list of samples used in the study is shown in Table 1. For both databases, the additional patient information related to demographic data, history of present illness, previous medical history, type of medication therapy, compliance with the therapy, as well as the level of hemoglobin A1c and other relevant laboratory values, was not available. Neither database contained information on the type of DM and if the normal looking images come from healthy individuals or from those diagnosed with diabetes that did not develop signs of DR.

| Automatic segmentation of the raw images
The raw images were segmented by using the software program developed as a part of the noninvasive medical imaging project focused on the diagnosis and early detection of non-communicable diseases. 28

| Selection of the region of interest
Diabetic retinopathy has a progressive course, and pathological changes close to macula, the area central vision, represent the most serious threat to eyesight. 5 Recent advances in OCTA allow visualization of the detailed 3-dimensional anatomy of vascular plexuses in healthy retina. OCTA showed that in the healthy retina the most dramatic changes in vessel density are present around foveolar avascular zone. 33 For these reasons, our research is primarily focused on the ROI microvascular network morphology in a circular area of retina that includes macula in the center and the optic disk at the periphery of each image. For the DIARETDB1 database, ROI was a circular area measuring 1000 pixels in diameter, containing the macula in the center and the complete optic disk. Since the images from the STARE database had a lower resolution and were centered inconsistently, the ROI was defined as a circular area of 500 pixels in diameter, centered either on the macula or on the optic disk. Each ROI was cropped using the GIMP image processor.

| Image binarization and filtering
Binarization was performed by using ImageJ software and its Mexican hat wavelet filter followed by adjusting the image threshold to maximum (white pixel). 34 The Mexican hat filter is commonly Box-counting dimension calculation is a type of monofractal analysis. The concept applied in multifractal analysis is analogous to applying filters to an image to emphasize certain features that otherwise are not obvious. 35 Structures that have monofractal nature are not affected by this process. While box-counting counts the number of boxes that contain foreground pixels, multifractal analysis counts the number of foreground pixels for each box and it assigns the value to the box according to that number, which is called the mass measurement. In multifractal analysis, ImageJ software sets an arbitrary exponent called Q to several values symmetrically distributed around 0 (eg, −3, −2, −1, 0, 1, 2, 3).
Following that, each mass measurement is emphasized by being raised to Q, and corresponding fractal dimensions denoted D −3 , D −2 , D −1 , D 0 , D 1 , D 2 , and D 3 are calculated. If the structure that is being analyzed has a multifractal nature, then D 0 >D 1 >D 2 . The f(α) curve is commonly used for interpretation of multifractal analysis. This is a convex curve, and its aperture from 1 to −1, as well as the aperture slope, can be used for overall assessment of the multifractal results that include for all values of Q. 35

| Statistical analysis
The t test and one-way ANOVA with Tukey post hoc test were performed by using the statistical program R. In order to measure the inter-rater variability, we calculated free-marginal kappa, a chanceadjusted measure of agreement for any number of cases, categories, or raters by using free online calculator. 36

| Automatic segmentation of images
The application of the automatic segmentation method produced gray scale images that specifically showed the vascular network details to high order of branching generations without artifacts from laser treatment, hemorrhages, and exudates. In addition, the use of two types of filters produced two sets of images. The application of the coarse filter resulted in a set of binarized vascular network images with a low number of vessel branching generations: vessel branches from the first to approximately the fourth generation of branching. The application of the fine filter yielded a set of binarized images with a high number of vessel branching generations that also included vessels belonging to more than four generations of branching ( Figure 3).

| Box-counting dimension and lacunarity of DIARETDB1 images with a high number of branching generations
Box-counting analysis of images with a high number of branching generations showed that the complexity of vascular network measured by the box-counting dimension Db increases as the severity of DR increases. At the same time, lacunarity (or gappiness, or inhomogeneity) decreases with more severe stages of DR ( Figure 4A,B).

| Multifractal dimensions of DIARETDB1 images with a high number of branching generations
The multifractal analysis of images with a high number of branching generations yielded results with the same trends as the results of the box-counting analysis ( Table 2): as the severity of DR increased, the complexity of the vascular network increased as well. In addition, this The analysis of f(α) spectra showed that the aperture length −1 to 1 is significantly shorter in the severe DR group, and the aperture slope has a more downward trend when compared to the moderate DR group ( Figure 5).

| Comparison of images with a high number and a low number of vessel branching generations: STARE and DIARETDB1 database
We The results of our previous report 27 support this observation. In that report, we used raw retinal images from the STARE database. It is important to emphasize that the resolution of these images was significantly lower compared to the DIARETDB1 images. Therefore, the level of details and the number branching generations on corresponding manually segmented images was relatively low. In that report, the analysis of these manually segmented STARE images showed that microvascular network in PDR displays increased lacunarity and decreased fractal dimension when compared to the healthy retinas. To explore this further, in the study we present here we performed the same automatic segmentation of the STARE images and applied the same filters as for the DIARETDB1 images. images with the high number of branching generations ( Table 3).
The multifractal analysis of images with low number of vessel branching generations did not demonstrate any significant changes in the microvascular complexity related to severity of DR (data not shown).

| D ISCUSS I ON
The analysis of DIARETDB1 images with a high number of branching generations showed that the increase in severity of DR is associated with increased complexity of microvascular network measured by box-counting dimension. At the same time, this is associated with decreased gappiness or inhomogeneity as measured by lacunarity.
In addition, we show that severe stages of DR are characterized by smaller aperture length of the f(α) spectra, when compared to moderate DR.
Taken together, our results demonstrate that microvascular  Automatic segmentation used in this study displays the ability to segment blood vessels without artifacts from hemorrhages, exudates, and laser photocoagulation. This allowed the fractal analy- of other systemic diseases, as well as diseases specific to the eye that could also be associated with changes in microvascular network morphology, was not available and is out of the scope of this study, but should be addressed in the future as well.

PER S PEC TIVE
Our results establish the fact that severe stages of non-proliferative and proliferative DR could be detected noninvasively by using a basic high resolution fundus photography and automatic segmentation of microvasculature to the high number of branching generations, followed by box-counting, lacunarity, and multifractal analysis of retinal microvasculature.
This approach shows that the increase in the severity of DR is associated with increased complexity of the microvascular network as measured by fractal dimension and decreased inhomogeneity of retinal microvascular network as measured by lacunarity.
Finally, we show that the results of fractal analysis are affected by the ability of a segmentation method to account for the smaller vessels with more branching generations.