Visual Tract Degradation in Bilateral Normal-Tension Glaucoma—Cortical Thickness Maps and Volumetric Study of Visual Pathway Areas

The aim of the study was to evaluate changes in the central visual pathways during the early and advanced stages of bilateral normal-tension glaucoma (NTG). Methods: The studied groups constituted patients with bilateral normal-tension glaucoma of the same stage (n = 45) and age-matched healthy volunteers (n = 17). All patients underwent ophthalmic examination and examination on a 1.5 Tesla Magnetic Resonance Scanner (Optima 360, GE Healthcare). Volume and cortical thickness analyses were performed using the open-source automated software package FreeSurfer. Results: There was a significant difference in lateral geniculate nuclei volume between the control and advanced glaucoma groups in the right hemisphere (p = 0.03) and in the left hemisphere between the early and advanced glaucoma patients (p = 0.026). The optic chiasm volume differed significantly between the control and advanced NTG groups (p = 0.0003) and between early and advanced glaucoma patients (p = 0.004). Mean cortical thickness analysis revealed a significant increase in values in the advanced glaucoma group in the right Brodmann area 17 (BA17) (p = 0.007) and right BA18 (p = 0.049) as compared to early NTG. In the left BA18 area, the mean thickness of the cortex in the early glaucoma group was significantly lower than in the control group (p = 0.03). Conclusions: The increase in the grey matter thickness in the V1 region with more-advanced glaucoma stages may reflect compensatory hypertrophy. Additionally, the regions of the brain early affected during glaucoma with reduced thickness were the right lateral occipital gyrus and left lingual gyrus. The most prominent change during the course of glaucoma was the increase in grey matter thickness in the right cuneus.


Introduction
Glaucoma is a heterogeneous group of diseases caused by the death of retinal ganglion cells (RGCs), and clinically characterized by the excavation of the optic nerve head and visual-field loss [1,2]. Glaucoma is the second-leading cause of blindness worldwide and the most frequent cause of irreversible blindness worldwide [2,3]. In Europe, the prevalence of glaucoma is 2.93% among persons aged 40 to 80 years. The prevalence rises with age, reaching 10% in persons over 90 years old [4].
The causes of RGC death in glaucoma remain unclear. Increased intraocular pressure (IOP) remains the most important and only modifiable risk factor for glaucoma development. However, there is an enormous group of patients who develop glaucoma, so-called normal-tension glaucoma (NTG), in the absence of IOP elevation. In the case of NTG, risk factors different from IOP elevations are suggested. Among them are a higher sensitivity to normal pressure, vascular dysregulation, an abnormally high translaminar pressure gradient, and a neurodegenerative process due to impaired cerebrospinal fluid dynamics in the optic nerve sheath compartment [5].
Generally, two principal theories for the pathogenesis of glaucomatous optic neuropathy have been described as a mechanical and a vascular theory [6]. Moreover, the studies suggest that glaucoma is similar to neurodegenerative diseases such as Alzheimer's disease and Parkinson's disease [7,8], with neurodegenerative changes localized in the central nervous system. The latter association is also supported by studies showing trans-synaptic degeneration affecting the central areas of the visual system in patients with glaucoma [7]. This shed new light on the disease, giving life to the new theory suggesting that the disease could be considered not only an ocular disease but a more complex neurodegenerative process that affects the entire visual system [7]. However, the involvement of the central nervous system during the course of glaucoma remains unclear, including whether the observed pathology is the cause or rather the result of the ongoing glaucomatous process.
The studies of the central involvement during glaucoma pathology also give another interesting possibility. Primary open-angle glaucoma (POAG) is painless and asymptomatic for a long time. Visible changes at the head of the optic nerve and defects in the visual field are detected when 40% of RGC is irreversibly lost [9]. For this reason, there is a clear need for the methods enabling the introduction of earlier detection. One of the promising techniques is magnetic resonance imaging (MRI) [10].
Normal-tension glaucoma (NTG), a subtype of POAG, occurs despite the IOP never exceeding the normal range. Recent studies have revealed that IOP-independent risk factors, including vascular factors, trans-laminar pressure gradient, immune-related disorders, genetic factors, or disturbances of the biomechanical properties of the eye, may play a crucial role in NTG pathogenesis [11,12]. Additionally, the relationship between glaucoma and cognitive impairment, especially that due to Alzheimer's disease and Parkinson's disease, was reported [13][14][15]. Most MRI studies did not separate NTG patients from the whole POAG group. However, patients with NTG were described to have abnormalities in the anterior visual pathway involving the optic nerve diameter, the height of the optic chiasm, and LGN volume [16]. Studies with the diffusion tensor imaging (DTI) technique indicated whole-brain white matter damage in NTG patients and showed a positive correlation with RNFL thickness [17,18].
In this study, we have applied both voxel-based morphometry (VBM) and surfacebased registration (SBR) methods for investigating anatomical alterations that occur in the brain during the glaucomatous process on a different degree of advancement. We measured volumes of the lateral geniculate nuclei (LGNs) and optic chiasm and cortical thickness in areas associated with vision using an open-source neuroimaging toolkit FreeSurfer (Massachusetts General Hospital, Harvard Medical School; http://surfer.nmr.mgh.harvard.edu, 27 January 2022).
The study aimed to evaluate the changes observed in the central visual pathways during the early and advanced stages of NTG.

Materials and Methods
The study was designed and conducted in the Department of Diagnosis and Microsurgery of Glaucoma and Department of Radiography in the Medical University of Lublin (Poland) with the acceptance of the Local Ethics Committee (approval number KE 0254/149/2018), and performed in accordance with the standards set by the Declaration of Helsinki. Before the study, all participants signed informed consent and underwent ophthalmic examination.

Participants
The studied groups contained patients with bilateral normal-tension glaucoma of the same stage (n = 45) and age-matched healthy volunteers (n = 17).
Patients were diagnosed with NTG according to the following criteria: glaucomatous neuroretinal rim loss, open angle in gonioscopy, and intraocular pressure (IOP) consistently <21 mm Hg (with the highest value ever measured in a given patient being ≤21 mm Hg, maximal IOP was evaluated by 4 measurements during office-hours), and no detectable eye pathology which would imply a diagnosis of secondary glaucoma. During the diagnostics, the maximal IOP was checked during 2 separate visits during different time points and, during every visit, it was measured twice with a 2 h break. Diagnosis was established by experienced ophthalmologists (EKJ, TZ).
The Humphrey visual field test 24-2 was performed on all patients The VF results were evaluated and the diagnosis was confirmed by glaucoma specialists (EKJ). Glaucomatous defect on SAP was defined based upon a glaucoma hemifield test result outside normal limits and the presence of at least 3 contiguous test points within the same hemifield on the pattern deviation plot at p < 1%, with at least 1 point at p < 0.5%, on at least 2 consecutive tests, with reliability indices better than 15%. According to the criteria of Hodapp-Parrish-Anderson in visual field examination, patients were divided into two subgroups: early glaucoma (MD better than −6 dB; n = 28) and advanced glaucoma (MD less than −12 dB; n = 17). Only the patients with bilaterally the same stage of the defect were included.
The control group was recruited from the patients who, after ophthalmic diagnostics, had excluded glaucoma and other diseases possibly influencing visual pathways (retinal or optic nerve disorders, central nervous system disorders). Demographic and clinical data of studied groups are presented in Table 1.

MRI Data Acquisition
All patients underwent examination on a 1.5 Tesla Magnetic Resonance Scanner (Optima 360, GE Healthcare) using a GE HD 8 channel neurovascular array coil. Standard brain and orbits protocol was extended by a 3D T1-weighted magnetization-prepared rapid acquisition gradient echo (MPRAGE) sequence acquired with the following parameters: TR = 9.4 ms, TE = 3.9 ms, FA = 8, TI = 1000 ms, FOV = 24 cm, matrix 192 × 192, slices = 166, voxel size 1.2 × 1.2 × 1.2 mm 3 . During 6 min acquisition, the patients were asked to remain still and to avoid any eye movement. MRI data acquisition was performed by experienced radiologists (RP, AP, PK).

Data Analysis
Volume and cortical thickness analyses were performed using the open-source automated software package FreeSurfer (version 7.1.1, Massachusetts General Hospital, Harvard Medical School; http://surfer.nmr.mgh.harvard.edu, 27 January 2022). Analyzed data were acquired using an automatic recon-all procedure in FreeSurfer containing the following steps: motion correction and intensity normalization, Talairach transformation, skull stripping and neck removing, registration, anatomical structure labeling, and white matter segmentation based on intensity, neighborhood, and smoothness constraints. Next, deformable surfaces models were created, starting from the tessellation process. Before generating final surfaces (white, pial, inflated) and registration to spherical atlas, automated topology fixer was used. Cortical thickness in that approach was calculated as the distance between white and pial surfaces.
Each subject underwent a quality assessment performed by an experienced radiologist. Prepossessing steps including brain mask and skull stripping check, normalization, and segmentation evaluation were executed manually to excluding the presence of topology defects and automated pipeline errors.
The volume of LGNs was measured in selected groups using a thalamic nuclei atlas developed by Iglesias et al. [19]. Optic chiasm volume assessment was performed according to data obtained from the Desikan and Killiany parcellation atlas [20].
Moreover, statistical cortical thickness maps were calculated using FreeSurfer's Query, Design, Estimate, Contrast (QDEC) graphical user interface. After smoothing cortical thickness measures using a full width at half maximum (FWHM) Gaussian kernel of 10 mm, each group and corresponding hemispheres were compared with each other using a General Linear Model (GLM) computed vertex-by-vertex. Results were presented with and without correction for multiple comparisons.

Statistical Analysis
PQStat software package ver. 1.8.2 was used for the statistical analysis of the mean volume and mean cortical thickness results. For group comparison, one-way analysis of variance (ANOVA) or Welch's ANOVA (in the case of unequal variances) was employed. For multiple comparison procedures, parametric Tukey's (when equal variances) and nonparametric Games-Howell's (when variances were unequal) post hoc tests were performed. When rejecting outliers, Chauvenet's criterion was applied. The results were considered significant at p < 0.05.
Significant results of cortical thickness change calculated using QDEC after clusterwise correction for multiple comparisons with a Monte Carlo Simulation (p-value threshold of <0.01) were presented on semi-inflated cortical surfaces. Results presented on maps were considered significant at clusterwise probability (CWP) < 0.05.

Results
Data obtained during thalamic parcellation and subcortical segmentation showed that there is a significant difference in lateral geniculate nuclei (LGNs) volume between the control (238.5 mm 3 ) and advanced glaucoma groups (208.2 mm 3 ) in the right hemisphere (p = 0.03), and in the left hemisphere between the early (240.5 mm 3 ) and advanced glaucoma (217.0 mm 3 ) patients (p = 0.026). Furthermore, optic chiasm volume differs significantly between the control (225.3 mm 3 ) and advanced glaucoma (156.9 mm 3 ) groups (p = 0.0003) and between early (206.2 mm 3 ) and advanced NTG patients (p = 0.004) ( Figure 1).
Mean cortical thickness analysis in anatomical regions selected from Brodmann's areas atlas revealed a significant increase in values in the advanced NTG group in the right BA17 (V1 area, p = 0.007) and right BA18 (V2 area, p = 0.049) as compared to early NTG. Moreover, in the left BA18 area, the mean thickness of the cortex in the early glaucoma group was significantly lower than in the control group (p = 0.03) ( Figure 2).
In the second analyzed anatomical atlas (Destrieux), the cortex was significantly thinned in the early stage of glaucoma as compared to the advanced stage in the cuneus (right, p = 0.013; bilateral, p = 0.034), left collateral sulcus and lingual sulcus (p = 0.034), and in the left occipital pole (p = 0.013). A significant increase in mean cortex thickness was observed in advanced NTG in comparison to the control group in the left occipital pole (p = 0.014). All results are presented in Table 2.
FreeSurfer's analysis performed in QDEC revealed visual regions from the Desikan and Killiany anatomical atlas of altered cortical thickness. Thinning of the cortex was observed in the early and advanced glaucoma groups in the right lateral occipital gyrus and in the left lingual gyrus as compared to the control group. Furthermore, there was an increase in cortical thickness values in the right pericalcarine of the advanced NTG group (compared to control). Evaluation of cortical thickness between both NTG groups showed an increase in values in the early stage in the right lateral occipital gyrus and left lateral occipital gyrus and a decrease in the right cuneus, right lingual gyrus, and right lateral occipital gyrus. Results are presented in Table 3.  Results after correction for multiple comparisons revealed significant thinning in the early glaucoma group as compared to advanced glaucoma only in the right cuneus (CWP = 0.045), in a cluster area of 431.1 mm 2 . Results are presented in cortical statistical maps in Figure 3.  Correlation analysis, presented in Table 4, revealed positive dependency between average temporal RNFL and mean LGN volume in a parametric test (r = 0.34, p = 0.027) in the early NTG group, and negative dependency between average nasal RNFL and mean LGN volume in a nonparametric test (R = −0.50, p = 0.035) in the advanced NTG group (compared bilaterally with corresponding eyes). When comparing opposite eyes, mean LGN volume correlated positively in the early glaucoma group with average RNFL    Table 5 shows the parametric (Pearson's) and nonparametric (Spearman's) correlations analyses based on clinical data (MD and average RNFL) and mean cortical thickness calculated in selected areas using FreeSurfer software. MD correlated negatively with V1 thickness in advanced NTG stage when comparing corresponding (r = −0.50, p = 0.009; R = −0.52, p = 0.00) as well as opposite sides (r = −0.40, p = 0.046; R = −0.48, p = 0.006). Furthermore, in the early NTG stage, correlations between MD and mean cortical thickness in the calcarine sulcus (R = 0.33, p = 0.034) and occipital pole (r = −0.30, p = 0.048; R = −0.35, p = 0.019) were noticed. There were no significant correlations observed between cortical thickness and average RNFL in any of the studied groups. Detailed results are presented in Table 5. Table 5. Parametric (Pearson's) and nonparametric (Spearman's) correlation analyses between mean cortical thickness in selected areas and the clinical parameters (MD and RNFL) in studied groups. Statistically significant results are marked bold.

Discussion
Glaucoma is a disease with complex pathogenesis in which the damage is detectable along the whole visual pathway. Degenerative changes were reported in the pregeniculate part, including the optic nerves and optic tracts, and in postgeniculate pathways with the optic radiations and visual cortex [24][25][26]. Additionally, recent studies have reported degeneration in non-vision-related areas of the brain, suggesting the presence of a visualpathway-independent brain involvement [27]. The two hypotheses explaining pathogenesis of the central involvement in glaucoma were proposed. One claims that the RGC loss leads to the transmission of the injury to the LGN and further to the visual cortex. However, there is also scant evidence showing that the LGN may be the primary site of the injury, with the secondary reduction in the axonal transport of neurotrophic factors causing RGC apoptosis [26,28]. The latter is proposed especially for NTG, where elevated IOP is not the main causative factor.
In its pathway to the brain, the optic nerve crosses the optic chiasm, so the LGN and primary visual cortex receive information from both the nasal visual field of the ipsilateral eye and the temporal visual field of the contralateral eye. This makes the assessment of the central visual pathway more difficult because the participation of the particular optic nerve in the central damage is unable to be estimated. To avoid this uncertainty, we decided to include only the patients with symmetrical binocular glaucomatous defects.
Changes in the optic chiasm are known in glaucoma. Zhang et al. examined the peripheral part of the visual pathway in NTG with MRI, and determined that constriction of the optic nerve and reduction in the size of the chiasm correlated with a decrease in RNFL [16]. Additionally, the decrease in the intracranial part of the optic nerve in glaucoma correlated with the histological decline in the number of RGC [29]. In this study, a significant decrease in optic chiasm volume was observed when glaucoma stages become more advanced. This difference reflects the proper inclusion of patients in the early and advanced glaucoma groups.
The LGN is the relay center in the thalamus for sensory information to enter the visual cortex. Furthermore, it plays a role in perception and cognition far beyond that of a relay nucleus, and it has to be considered as an early gatekeeper in the control of visual attention and awareness [30]. The LGN consists of six distinctive layers. The inner two layers are magnocellular layers, while the outer four layers are parvocellular layers. In between these two layers are koniocellular cells. Numerous studies using MRI imaging show changes in LGN in the course of glaucoma [31][32][33][34]. Similar changes were observed in animal models of glaucoma [35]. Our study in NTG patients showed that there is a significant difference in LGN volume between the control and advanced glaucoma group in the right and in the left side between the early and advanced glaucoma patients. Additionally, in early glaucoma the decrease in LGN volume was correlated with RNFL thinning measured with OCT. Similar to our previous results [25], such correlation was not observed for advanced glaucoma. Although the general tendency to decrease in LGN volume in more-advanced stages of glaucoma is observed, the relationships seem not to be very distinctive. There might be some reasons for this phenomenon. First, two main kinds of neurons in LGNs are described: interneurons confined to the LGN, and relay neurons with synapses to the visual cortex. Interneurons were suggested to be more resistant to transneuronal degeneration [32,36,37] during the glaucoma process. Similarly, two main types of pathways in LGN, magno-and parvocellular, are claimed to have different involvement during glaucoma. Moreover, the early glaucoma group may be the most heterogeneous, since before the VF changes occur the prominent RGC percentage needs to be injured. On the other side, our results may reflect the neurodegenerative process occurring independently in LGNs during NTG. The studies regarding neurodegenerative changes in central structures in the course of glaucoma are mainly conducted on animal models with high-tension glaucoma (HTG) and do not provide knowledge of the processes operating during NTG in humans. Advanced MRI techniques may help to obtain a clearer view regarding central involvement in human NTG.
The human cortex is a highly folded structure, so measuring cortical thickness based on MRI data requires complex reconstruction algorithms. The average thickness of the cerebral cortex for the whole brain is approximately 2.5 mm [38]. Obtaining MR images with submillimeter resolution is time-consuming and not always possible in clinical conditions. Typically, 3D T 1 -weighted images used for measurements are characterized by a resolution of 1.5 mm 3 or less. To address these problems in our study, we used a technique called the surface-based registration (SBR) method that defines cortical thickness as the distance between two surfaces: grey/white and pial [39]. It allowed for precise calculation independent from voxel size and obtained results with accuracy around 10 −4 m [39].
While analyzing the thickness of the cerebral cortex, we used three anatomical atlases [20][21][22][23], each one of them defining the areas of the visual cortex somewhat differently. The Destrieux atlas is the most detailed as compared to the Brodmann areas and Desikan and Killiany atlases. For example, only the Destrieux atlas contains an area called the occipital pole, while in the Desikan and Killiany atlas this region is included as a small part of a few neighboring labels. We aimed to check cortical abnormalities in NTG glaucoma in the most exact way possible, so we have decided to evaluate results according to different regional patterns.
The primary visual cortex is the predominant brain area receiving direct input related to visual stimulation. V1 is the first of the cortical regions to process information and is divided into six layers, each with different cell-types and functions. V1 responds to simple visual components such as orientation and direction [40]. V2 receives integrated information from V1, sends feedback connections to V1, and has feedforward connections with V3-V5. In our study, cortical thicknesses of V1 and V2 differed between the studied groups. The decrease (not significant) was observed in early glaucoma compared to control and, surprisingly, the increase in the cortical thickness was described in the advanced group. Similar increased volume of the components of the visual system, which enables for identification and classification of visual stimuli, was described by Williams [41]. The lack of statistical difference in V1 thickness between early glaucoma and control groups was reported previously for POAG [42], and was explained by the fact that the majority of retinal inputs project to layer 4 of the V1 area [43], receiving information from LGN. Hence, the possible selective degeneration of layer 4 in the V1 area may not significantly influence the whole V1 thickness [42].
V1 sends axons to the visual association areas [44]. Therefore, it could be speculated that degeneration of the visual pathway [45,46] could lead this region to increase its functional segregation ability in a compensatory capacity. According to our results, these compensatory mechanisms are observed in advanced, but not early glaucoma, which will support the hypothesis that the hypertrophy is the compensatory mechanism for visual-pathway atrophy rather than that temporary compensatory hypertrophy is followed by final atrophy. Our patients from the advanced group had bilaterally really advanced glaucoma with MD values near legal blindness. Increased cortical volume may be a sign of neuronal injury with cellular swelling or microglial activation [47,48]. Conversely, it may also reflect enhanced neuronal function or cortical plasticity [49,50]. There is growing evidence that extensive reorganization in response to sensory deprivation is possible in the adult brain and the anatomical changes may result from these functional adaptations. In the case of blindness, many studies have demonstrated that the visual association cortices are recruited for diverse functions including tactile, auditory, and linguistic tasks [51,52]. In this study, the tendencies similar to those observed for V1 were present also for cuneus and precuneus, areas classically related to visual information processing [53].
Fusiform gyrus and lateral occipital gyrus are located in the visual association cortex (V2), whose major function is the integration of visual information and generation of conscious perceptions [54]. It is said that V2 degenerates rapidly while V1 changes little during glaucoma progression. The reason is that V2 neurons may possibly be modulated by more-complex properties than V1, such as visual form [55,56] and motion processing [57,58], both of which are affected during POAG [42]. Additionally, V2 neurons are modulated by complex properties such as visual form and visual illusion [59], which may be impaired by the progression of the disease. Wang et al. [54] described that in HTG patients the V2 degenerated significantly while the primary visual cortex was relatively untouched. In NTG, we observed a significant increase in lateral occipital pole in advanced glaucoma in the left side, whereas there was no change in the right side. Similar lateralization to the left side of the damage was described previously [42,[59][60][61][62], however, in our results some changes were more prominent for the right side. Our patients had a similar degree of glaucomatous damage in both eyes. Slightly different functions for each of the brain hemispheres were described, for example, face recognition [63], one of the tasks affected by glaucoma [64][65][66]. On the other hand, further studies on patients with similar stages of glaucoma bilaterally evaluating the influence of lateralization are needed to exclude bias.
To confirm the regions of the brain affected by glaucoma, FreeSurfer's analysis was performed in QDEC, which revealed visual regions from the Desikan and Killiany anatomical atlas of altered cortical thickness. When compared to the control, the regions affected in early glaucoma with reduced thickness were the right lateral occipital gyrus and left lingual gyrus. These regions were also affected further during the course of glaucoma. The correlation between thinning of visual cortex and clinical measurements was previously described for glaucoma [42,59,61,67]. Correlation of gray matter volume with glaucoma severity was also demonstrated by significant differences in visual cortex thickness between mild and severe glaucoma groups, with increased damage in late-stage glaucoma [59,68]. In our study, when comparing early and advanced groups the most prominent difference was an increase in the GM thickness in the right cuneus during the course of glaucoma. One of the cuneus functions seems to be to integrate the somatosensory information with other sensory stimuli and cognitive processes such as attention, learning, and memory [53]. These results indicate that NTG predominantly involves visual cortices, accompanied by abnormalities of other cortices associated with attention and cognition. However, the interactions between these cortices remain unclear. NTG changes the information processing in the brain by reorganizing the information flows between visual cortices and non-visual brain areas, and these alterations are consistent with microstructural injury of visual cortices [69]. Li reported decreased information outflows in NTG patients from the left BA17 to bilateral cuneus and increased effective connectivity, reflecting integration within a distributed system, from BA19 and BA18 to cuneus (BA17) [70].
When evaluating patients with glaucoma, the structure (RNFL) and function (VF) of the visual pathway are measured and, usually, they are not strictly related. The assessment of early and late stages of the disease is particularly challenging with only structural defects in early phase and functional progression accompanied by stable RNFL thickness in advanced stages. The interesting finding in our study is the negative correlation between MD value resembling the depth of functional deficit and the thickness of BA17 observed in advanced glaucoma, which may support the hypothesis that in advanced glaucoma the measurement of visual cortex thickness may be an additional diagnostic tool to assess possible progression. Qing et al. [71] also showed that POAG neuropathy leads to decreased cortical activity in the primary visual cortex corresponding to the central normal visual field, but not OCT measurements and ophthalmoscopy. In no glaucoma stage were the visual cortices measurements correlated with RNFL thickness. Such correlations were present only for early glaucoma and LGN volume, the first cortical visual center.

Conclusions
To sum up, in this study we evaluated changes in the brain regions connected with the visual functions in different NTG stages. We showed the increase in the GM thickness in the V1 region with more-advanced glaucoma stages correlated with VF parameters, which may show compensatory hypertrophy. Additionally, the regions of the brain affected early during glaucoma with reduced thickness were the right lateral occipital gyrus and left lingual gyrus. When comparing early and advanced groups, the most prominent difference was increase in the GM thickness in the right cuneus during the course of glaucoma.
Author Contributions: T.Ż., R.P. and E.K.-J. are responsible for conceptualization; methodology was designed by A.P., S.M., P.K., T.Ż., R.P. and E.K.-J.; initial investigation was conducted by A.P., S.M., P.K. and E.K.-J.; A.P. and P.K. performed experiments; data were analyzed by A.P.; data visualization was conducted by E.K.-J. and A.P.; the original draft was written by: A.P., S.M. and E.K.-J.; resources were provided by E.K.-J., S.M. and T.Ż.; writing-review and editing and supervision were performed by T.Ż. and R.P. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.
Institutional Review Board Statement: The study was designed and conducted with the acceptance of the Local Ethics Committee (approval number KE 0254/149/2018) and performed following the standards set by the Declaration of Helsinki.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

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