Local neuroplasticity in adult glaucomatous visual cortex

The degree to which the adult human visual cortex retains the ability to functionally adapt to damage at the level of the eye remains ill-understood. Previous studies on cortical neuroplasticity primarily focused on the consequences of foveal visual field defects (VFD), yet these findings may not generalize to peripheral defects such as occur in glaucoma. Moreover, recent findings on neuroplasticity are often based on population receptive field (pRF) mapping, but interpreting these results is complicated in the absence of appropriate control conditions. Here, we used fMRI-based neural modeling to assess putative changes in pRFs associated with glaucomatous VFD. We compared the fMRI-signals and pRF in glaucoma participants to those of controls with case-matched simulated VFD. We found that the amplitude of the fMRI-signal is reduced in glaucoma compared to control participants and correlated with disease severity. Furthermore, while coarse retinotopic structure is maintained in all participants with glaucoma, we observed local pRF shifts and enlargements in early visual areas, relative to control participants. These differences suggest that the adult brain retains some degree of local neuroplasticity. This finding has translational relevance, as it is consistent with VFD masking, which prevents glaucoma patients from noticing their VFD and seeking timely treatment.

www.nature.com/scientificreports/ However, whether the mere presence of differences in size and position of RFs is evidence of cortical reorganization has been questioned, as similar changes occur with simulated (artificial) scotomas 15,17,18,[26][27][28][29] . Consequently, it remains unresolved whether such changes in representation are the result of reorganization or simply the consequence of the damage at the level of the eye changing the input that reaches the cortex. Therefore, to fully address the issue of retained adult plasticity in glaucoma, it is essential to include studies on participants with peripheral VFDs, to use high-precision analyses, and to include control conditions that account for the visual deprivation due to natural scotomas.
For these reasons, in this study, we investigated how glaucoma affects the functional organization of the visual cortex using fMRI in combination with advanced neural modeling. We assess the fMRI responses of participants with glaucoma vis-a-vis those of age-matched control participants with a matched, simulated VFD. This is crucial, as it will allow us to rule out that any observed changes in the pRFs are merely the result of the altered visual input 15,17,25,30 . Therefore, for each participant with glaucoma, a matched control participant observed the visual stimuli with a simulated scotoma (SS) designed to mimic the glaucoma participant's reduced visual sensitivity, as assessed using Standard Automated Perimetry (SAP). Moreover, we assessed retinal thickness using Optical Coherence Tomography (OCT) and applied fMRI-based VF mapping techniques, based on both standard pRF mapping and an advanced variant called micro-probing 31 . The latter approach better predicts the pRF shape and is more accurate when charting the VF via fMRI 32,33 .
To preview our results, while the coarse retinotopic reorganization of the visual cortex was maintained in the participants with glaucoma, we found a reduction in its Blood Oxygen Level Dependent (BOLD) signal responsiveness and differences in the distributions of estimated pRF sizes and positions compared to those in control participants with a matched SS. In other words, the observed pattern of reorganization of the visual cortex in glaucoma appears not to be a mere consequence of the reduced visual input due to the retinal scotoma. Moreover, when expressed in terms of the changes in pRF size and location, we find that the degree of reorganization correlates with the severity of the glaucomatous VF damage. In addition, fMRI-based VF reconstructions showed that glaucoma participants exhibit a lower VF sensitivity compared with controls with SS and local differences between fMRI-based VF sensitivity patterns from those obtained via SAP. Together, our findings suggest that the adult visual cortex retains a spatially localized capacity to functionally reorganize.

Results
Glaucoma affects the BOLD modulation in the visual cortex. Figure 1A shows the BOLD modulation as a function of eccentricity for participants with glaucoma and for the control participants with and without a simulated scotoma, controls SS and controls NS, respectively. At all eccentricities, the BOLD modulation is reduced in the participants with glaucoma compared to both control conditions (F(2,31) = 12.98, p < 0.0001, FDR-corrected). There was no significant evidence for differences in the BOLD modulation between visual areas, V1, V2 and V3 (F(2, 31) = 2.67, p = 0.07) nor for an overall effect of eccentricity (F(12, 341) = 0.49, p = 0.92). The slightly higher foveal responsiveness of the control participants in the SS compared to the NS condition caused a significant interaction of group and eccentricity (F(12, 341) = 2.05, p = 0.02). Figure 1B,C show the average BOLD modulation as a function of macular nerve fiber layer (MNFL) thickness and contrast sensitivity, respectively. Data are plotted per individual VF quadrant. For both parameters, we find significant correlations with average BOLD modulation (contrast sensitivity: r 2 = 0.31, p = 0.006 ; MNFL thickness r 2 = 0.31p = 0.007 ). As shown in Fig. S1, the MNFL thickness correlation is also present for visual areas V2 ( r 2 = 0.34p = 0.0024 ) and V3 ( r 2 = 0.35p = 0.0017 ) while a significant correlation of BOLD vs contrast sensitivity is present for V3 but not for V2 (V2: r 2 = 0.2, p = 0.08 ; V3: r 2 = 0.22, p = 0.05 ). Figure 1D shows the V1 BOLD modulation of control participants in the SS condition as a function of simulated contrast sensitivity. The correlation is not significant ( r 2 = 0.08, p = 0.44 ). As shown in Fig. S1, this was neither the case for V2 ( r 2 = −0.01, p = 0.92 ) nor for V3 ( r 2 = 0.09, p = 0.42). Figure 2A shows the pRF properties (eccentricity, polar angle and size) projected on the inflated brain mesh, obtained for participant pair P03 consisting of a participant with glaucoma (G03) and matched control participant (C03). The latter performed the experiments both without (NS) and with (SS) a simulated scotoma matched to the scotoma of G03. Despite the severe glaucomatous VFD, represented by the VF plots measured using SAP in panel 2B, the global retinotopic organization of the visual cortex is preserved in G03. Figure S2 shows the retinotopic maps of additional glaucoma participants with different VFDs. Panels 2E and 2F show that the eccentricity and polar angle histograms exhibit the same pattern for G03 and C03 in both the NS and SS conditions, and only local differences can be observed. Regarding pRF size, overall we found no significant differences between G03 and C03 SS and NS. Still, as shown in Panel 2G, at more peripheral locations (5° to 7°) which are more affected by the VFD both G03 and C03 (SS condition) show larger pRFs sizes than C03 (NS condition). As expected, G03 and C03 (SS and NS) show an increase of pRF size with eccentricity. Panel 2D shows various single-voxel time series for the glaucoma and control participant. The voxel locations are indicated by the colored arrows in panel 2A. The top graph of panel 2D shows the time series recorded in two voxels. One voxel is located in the Scotoma Projection Zone (SPZ) of G03 (red arrow) which, in this case, is situated in the periphery of the lower left quarter field. Of relevance, according to SAP, G03 still had residual sensitivity within their SPZ. The second time series was recorded in a voxel with a (mirrored) pRF position in the upper left visual quarter field. Therefore, this was located in the contralateral hemisphere (orange arrow). Not so surprising, the BOLD modulation of the voxel inside the SPZ (dark red line) is substantially lower than the one located in the contralateral hemisphere (orange line). Moreover, this decrease was substantially larger than expected based on the reduced stimulation that is the consequence of the VFD, as mimicked by the simulation in C03 (SS), as shown in the bottom graph of panel Figure 1. V1 BOLD modulation varies between participants with glaucoma and control participants and correlates with disease severity. (A) V1 BOLD modulation, defined as a function of eccentricity for glaucoma (red) and control participants with (dark blue) and without a simulated scotoma (light blue). The BOLD modulation was binned in 1° eccentricity bins. The bars represent the 95% confidence interval (CI). (B) Correlation of the V1 BOLD modulation of participants with glaucoma obtained from individual quadrants with MNFL thickness (calculated by averaging the values of the macula of both eyes). Each data point is from a separate quadrant of an individual participant with glaucoma. The different shades of red and blue colors denote the different Glaucoma and Control participants, respectively. (C,D) Correlation of the V1 BOLD modulation from separate quadrants with the contrast sensitivity of both eyes combined (the max between the mean deviation of the two eyes) for glaucoma and controls, respectively. Each data point is from an individual VF quadrant. The vertical dashed lines and shaded blue area represent the average values ± SD of the MNFL thickness (plot B) and contrast sensitivity (plots C,D) measured in control participants, Table 1. The correlations between the BOLD modulation and the disease severity (contrast sensitivity and MNFL thickness) were calculated using a linear mixed effects model with a slope and intercept per subject as a random effect. Detailed information regarding the correlation analysis is provided in "Correlation analysis". www.nature.com/scientificreports/ pRF properties in the inflated brain mesh for G01 can be found in Fig. S2. Panels 3A and 3B indicate that while for C01, the pRF distributions in the NS and SS conditions are similar, for G01 the normalized number of pRFs located in the scotomatic region (the binocular scotoma of G01 is primarily located in the upper left quadrant defined by a polar angle between 0 and π 2 , Panel 3D) is reduced compared to that in the SS condition of C01. Figure S3 shows the normalized distributions for the remaining participant scotoma-control pairs. Furthermore, Fig. 3C shows that G01 exhibits larger pRFs than C01 which are not restricted to the scotomatic area but are present throughout the entire (central) VF. For completeness, we present the analysis of the variation of the pRF size as function of the polar angle in Fig. S4, which also shows that the glaucoma participant has larger pRFs sizes across the entire visual field when compared to the control participant, but with and without the simulation (SS and NS). However, such larger pRFs were not consistently found in all glaucoma participants. For example, in G03 depicted in Fig. 2, the pRF size maps provide little evidence for deviations in pRF size.

Large scale organization of the visual cortex is preserved in glaucoma.
In fact, as we will report further on, heterogeneity rules. PRF size can vary substantially even in healthy controls, for example due to eccentricity, eye movements, differences in SNR and presence of simulated scotomas. We take this into account in our analysis by comparing pRF size of Glaucoma vis a vis to the baseline deviation of the control SS versus all of the controls NS. Thus the heterogeneity in pRF size measure in glaucoma participants is larger than in controls. We compared the normalized distribution of the pRF estimates between each glaucoma participant and their respective control. While we found deviations in all glaucoma participants, in some this primarily affected pRF size, while in others this was mainly reflected in the pRF location (see Figs. S7, S8 and S9). In particular Fig. S7, shows that glaucoma heterogeneously affects pRF size, evident from multiple observable patterns. The majority of glaucoma participants (10 out of 19) show larger pRFs compared to the matched controls (evident from a negative deviation for smaller pRFs and a positive deviation for larger pRFs, e.g. P01). Others show the opposite pattern and have smaller pRFs compared to the matched controls (evident from a positive deviation for smaller pRFs and a negative deviation for the larger ones, e.g. P19). Yet, others show no significant differences at all (e.g. P15). However, when grouped together, on average, there is no significant difference between the V1 pRF size of participants with glaucoma and the control participants with and without SS (t(60)) = 1.63, p = 0.10; Fig. 3E). Figure S3 shows that also for visual areas V2 and V3 there were no such differences (t(60)) = − 0.65, p = 0.52 and t(60)) = − 0.26, p = 0.80, respectively).
In order to perform a group analysis and taking into account the heterogeneity of the VFDs of glaucoma participants, we correlated the average pRF size per quarter field with the metrics obtained via ophthalmic tests (contrast sensitivity measured with SAP and MNFL thickness measured with OCT). Panel 3F shows that for the participants with glaucoma, the average V1 pRF size per quarterfield did not correlate with macular thickness (r 2 = 0.008, p = 0.7) while panel 3G shows, again for the participants with glaucoma, that the average V1 pRF size per quarterfield did correlate with SAP VF sensitivity (r 2 = − 0.28, p = 0.02), with larger pRFs in case of lower contrast sensitivity. Figure 3H shows that in controls with simulated visual field defects, this effect was not found (r 2 = − 0.02 p = 0.84). Similar patterns are present for visual areas V2 and V3 (Fig. S4). We note that for the glaucoma participants, the outcomes were substantially influenced by the three data points of three different participants (G01, G03 and G04) for which contrast sensitivity is below − 10 dB. Removing these from the analysis leaves the correlation between pRF size and SAP significant for V3 (r 2 = − 0.25; p = 0.03) but not for V1 and V2 (V1: r 2 = − 0.05; p = 0.7; V2: r 2 = 0.09; p = 0.43). These results suggest pRF enlargement might be a possible mechanism to compensate for the effect of a (severe) VFD, which would be in line with the finding of increased spatial integration as found with psychophysics 34,35 .
Variation in average pRF location between Glaucoma participants and controls SS is likely associated with contrast sensitivity loss. Figures S8 and S9 show the deviations in pRF eccentricity and polar angle between participants with glaucoma and the respective control participants. In all 19 glaucoma-control (SS) participant pairs, the distributions differed significantly from baseline in at least one of the bins. However, the location of these differences was highly variable between participants. Consequently, when evaluated at the group level, these differences tend to cancel each other out. This is confirmed in Fig. 4A,B, which show that when aggregated across all the participants with glaucoma and the control participants (SS), there are hardly any differences in the average normalized histograms of pRF eccentricity and polar angle. Importantly, at the individual level, we show that Euclidean Distance (ED) between the average position calculated per quarter field for glaucoma and respective control SS participants correlated significantly with contrast sensitivity (r 2 = − 0.23; p = 0.049, panel D) but not with macular thickness (r 2 = 0.08; p = 0.48, panel C). The within-participant ED between the SS and NS conditions does not correlate significantly with the simulated loss of contrast sensitivity (r 2 = − 0.00003; p = 0.99, panel E). Similarly to the relation between the pRF size and contrast sensitivity, also the relation between ED of glaucoma and Controls SS is weighed by the three data points with the lowest contrast sensitivity. Without these three points, the correlation becomes non-significant (r 2 = 0.022; p = 0.85). Furthermore, this analysis is corroborated with a more complex, detailed and locally specific analysis, based on local deviations of the pRF properties distributions between participants with glaucoma and control participants SS. In line with Fig. 4, Fig. S10 shows that the pRF position deviations between participants with glaucoma and control participants SS correlated significantly with contrast sensitivity (r 2 = − 0.24; p = 0.004, panel A) but not with macular thickness (r 2 = − 0.04; p = 0.74, panel C). The deviation between control SS and NS does not correlate significantly with the simulated loss of contrast sensitivity (r 2 = − 0.1; p = 0.39, panel B). Figure S3 shows that the pRF parameter distributions are consistent across participants.
Glaucoma affects the visual system beyond the eye. To investigate whether: (a) the brain representation of the VF matches the one measured at the level of the eye and (b) there is a higher sampling of the SPZ that could explain some perceptual phenomena characteristic of glaucoma, such as predictive masking, we Nevertheless, this reduction is not as strong as expected based on SAP, in particular in the periphery. In addition, although to a lesser extent, also the control SS participant shows a reduction in VF sensitivity in the quadrants most affected by the SS. Importantly, such reduction in VF sensitivity is not present in the control participant (NS). Furthermore, the VF reconstruction results are also in line with the pRF analysis shown in Fig. 2. This confirms that: (1) the VF reconstruction techniques are accurate and allow to retrieve VF sensitivity in glaucoma; and (2) that the SS have the desired effect although at a smaller scale. Importantly, the real and SS scotomas become smaller with visual hierarchy which can explain the predictive masking of the scotomas experienced by glaucoma patients. Note that the, SS can also introduce some biases in the VF reconstructions, in particular in visual areas whose SNR is low, i.e., the VF reconstructions of the Controls SS for V2 and V3 contain some reduced sensitivity regions in the lower VF where no reduction in contrast sensitivity was simulated. When examining the VF reconstructions ( Fig. 6) across the pairs of the participants with glaucoma and their controls with matched simulation (SS) it is clear that the VF sensitivity in glaucoma participants is overall lower than in controls SS. The majority of the controls SS do not show a marked scotoma. This is related to the fact that: (a) within the FOV of stimulation within the scanner only one glaucoma participant has a binocular scotoma with a sensitivity below − 15 dB (which is the threshold for being considered blind in that portion of the VF), and (b) a reduction in contrast does not significantly affect the pRF estimates 36 . The differences between the glaucoma and control SS participants suggest that either SS based on contrast sensitivity measured by SAP cannot accurately reproduce the real scotomas, or that the damage caused by glaucoma goes beyond the eye and affects the brain. Secondly, different patterns are observed, while some glaucoma participants (i.e., P2, P5, P6, P7, P16, P17) exhibit a lower VF sensitivity mapped via fMRI compared to what was measured with SAP, while for others the reverse pattern in observed (i.e., P1, P3, P4, P8, P11, P19). This suggests that glaucoma participants develop different strategies to cope with their VFD and that the effect of glaucoma goes beyond the eye. Some controls SS show increased VF sensitivity compared to the baseline controls NS (green regions in the VF reconstruction graphs), these increase VF sensitivity may reflect first the variation in fMRI based VF reconstructions and the effect of SS inducing short-term pRFs shifts directed to perceptually mask the scotoma 37 .

Discussion
In this study, we used fMRI in combination with model-driven analyses to quantify changes in cortical functional organization in participants with the ophthalmic disease glaucoma. Our findings suggest that visual cortex functioning is altered in glaucoma, which could be indicative of cortical plasticity. The main findings were: (1) cortical areas retained their coarse retinotopic organization, consistent with what has been described previously 15,38 ; (2) marked differences in the magnitude of the BOLD response that cannot simply be explained by reduced contrast sensitivity, (3) local differences in the size and position of pRFs, and (4) notable differences in the fMRI-based VF reconstructions obtained for participants with glaucoma compared to those for matched control participants with similar but simulated visual loss. Our findings suggest that glaucoma is associated with limited and local functional remapping in the early visual cortex. These local neural reconfiguration patterns are highly variable between participants with glaucoma and likely depend on the preserved VF. Such local changes may be interpreted as attempts by the visual system to compensate for loss that results from the neural damage. Moreover, the observed pRF changes are consistent with predictive masking of the scotomas and may have implications for future treatments. Below, we discuss our findings and interpretation in detail.
In participants with glaucoma, throughout their entire VF, the magnitude of the BOLD modulation is reduced when compared to control participants with matched simulated VFDs. The comparison to the latter indicates that these changes cannot be explained by reductions in contrast sensitivity reducing the strength of the signals reaching the cortex. Importantly, the reduced BOLD modulation in glaucoma participants takes place despite the higher number of retinotopic mapping runs acquired for glaucoma participants (6 runs) compared to control participants with (SS, 4 runs) and without a simulation (NS, 4 runs). This larger number of retinotopic mapping runs for glaucoma participants was included because we anticipated that this would be required to obtain reliable pRFs with deteriorated vision in some parts of the VF, and a difficulty to fixate for some of them. Moreover, when analyzed per quadrant of the VF, the modulation of the fMRI signal correlates with both the contrast sensitivity and the macular thickness per quadrant. In contrast, in controls with simulated VFDs this correlation is absent. Moreover, our results support the occurrence of glaucomatous alterations beyond the scotomatic regions, apparent from decreased BOLD activity throughout the primary visual cortex. These functional alterations extend beyond V1 and affect also V2 and V3. Our findings may be explained by transsynaptic degeneration, in which damaged nerve fibers in the optic nerve, and the subsequent death of retinal ganglion cells, will affect the entire visual pathway (including the optic nerve, lateral geniculate nucleus, optic radiation) and visual cortex, resulting in fairly widespread neuronal loss [39][40][41] . Our finding of reduced BOLD modulation in glaucoma participants agrees with previous studies [42][43][44] . Only one previous study found no association between the BOLD signal and RNFL thickness 45 . The reduced BOLD amplitude in early visual cortex may also be the result of deficient cortical perfusion. This explanation would be in line with previous studies that showed that participants with glaucoma, when compared to controls, have reduced cerebral blood flow in early visual cortical areas 46,47 . Analogous to this issue at the level of the retina, the question remains whether such limitations in perfusion cause reduced functionality and evoke neural degeneration, or merely reflect reduced oxygen consumption by degenerated neural tissues 48  www.nature.com/scientificreports/ Based on a retinotopic mapping analysis we found that, overall, in the participants with glaucoma, their cortex retains its coarse retinotopic organization. Locally, however, their cortices have deviant neural configurations as evident from differences in sizes and positions of the pRFs. Importantly, these differences were present in comparison to control participants with matched simulated scotomas. Yet, the pattern of changes varied substantially between participants with glaucoma. Given that our glaucoma population has heterogenous VFD that are located at different positions within the VF and are of variable extent, a group level analysis would be insensitive to any local scotoma specific changes. This, therefore, called for a participant-specific analysis, in which the differences in pRF estimates within each glaucoma-control pair were compared to the expected baseline intersubject variability in pRF properties observed in the controls. This analysis, presented in detail in the SI, suggests that each glaucoma participant is unique in their local pRF reconfigurations. While for some participants the reconfiguration manifests in changes in pRF size, in others it manifests in different pRF positions.
Overall, the size of pRFs and the deviations in their position distributions tend to increase with larger damage to the VF (i.e. with a larger loss in contrast sensitivity). Notably, we did not observe the same effect in the controls with simulated scotomas. We do note that the statistical significance of the relation between pRF size and positional deviations on the one hand, and the severity of the VF damage on the other hand, depends on the inclusion of the few data points with lowest contrast sensitivity (< − 10 dB) of different participants. To corroborate this relationship, it would be necessary to acquire data from more participants with severe VF damage in binocular quarter fields. This is important given that the pRFs in the real and simulated scotomatic regions may be modeled on signals with a lower SNR. Therefore, these could be more affected by modeling biases that result in artifactual pRFs enlargement 36 . To better assess this, we calculated the relation between the variance explained of the pRF estimates and the contrast sensitivity. Despite an overall lower SNR of the glaucomatous BOLD signal, which likely underlies less accurate pRFs (lower VE, Fig. S16), there is no bias associated with contrast sensitivity. Of relevance, the pRFs models estimated in severely damaged quarters of the VF present a VE that is in the range that is measured in the healthy VF (0.4-0.6, see Fig. S16). Moreover, also in the controls with simulated scotomas, we found no significant association between the VE of the pRFs models and the simulated contrast sensitivity. Nonetheless, in controls, a comparison of the pRFs estimated in the presence or absence of severe simulated scotomas, indicates that the former have a lower VE. These are thus more prone to reflect analytical artifacts which can be misinterpreted as short-term plasticity. In future work, it would be highly worthwhile to perform an extensive simulation study to determine the influence of SNR in interaction with SS on the estimated pRF properties, including evaluating multiple SS contrasts and sizes. Despite these factors, we interpret the differences in pRFs properties as evidence for limited local cortical plasticity in adults with glaucoma. These local changes in pRF properties are consistent with what was previously described in homonymous VFDs 15,38 . This reorganization may involve the activation of long-range horizontal connections in the visual cortex, such that healthy neurons in the cortex surrounding the lesion thus take control of the deprived ones. These reconfigurations enable the neurons within the lesion projection zone to capture information from spared portions of the VF 4,11,12,14 . This not only holds true for the visual cortex, but also for other sensory areas of the adult cortex [49][50][51][52] , cortical maps can for instance also reorganize after loss of sensory afferent nerves from a limb 53 . Indeed, such capturing of information from outside the lesion projection zone would be required for predictive masking of the natural scotomas to occur, a phenomenon that is clinically frequently observed 22 .
The VF reconstructions based on the fMRI data show that the cortical sensitivity of the glaucoma participants is reduced compared to the controls with simulated scotomas. This suggests that the glaucomatous deterioration of the visual system goes beyond the retina or that the way that the visual input is processed at the level of the retina is altered in glaucoma. In addition, while some glaucoma participants exhibit a lower VF sensitivity mapped via fMRI compared to what was measured with SAP (including damage in regions of the VF that appear to be 100% functional in SAP) others seem to show much smaller scotomas than what was measured via SAP. There are multiple explanatory mechanisms for these differences: (1) a reconfiguration of the RFs may cause neurons that were initially located inside the scotoma projection zone to process information from the spared VF, resulting in predictive masking of natural scotomas, and therefore the scotomas become smaller than measured in via SAP 54 , (2) fMRI allows to detect subtle changes in the visual system functioning to which SAP is insensitive and (3) the accuracy of standard automated perimetry (SAP) measures was insufficient. Furthermore, the scotoma projection zone shrinks across the visual hierarchy. This finding supports the view that predictive masking of the scotoma (the reason why glaucoma participants cannot perceive their scotomas), may result from feedback from higher cortical areas where scotoma representation is small or inexistent to earlier visual areas. This mechanism can also be the driving force behind pRF shifts 37 .
Various of our findings suggest that visual cortex functioning is altered in glaucoma. Glaucoma participants show reduced BOLD modulation compared controls; pRF size and position deviations tend to be larger in quarterfield sections that showed greater loss in contrast sensitivity, and VF sensitivity of glaucoma participants differ from controls SS. These cerebral adaptations have important applications as the clinical diagnosis and treatments for glaucoma are currently only focused on the eye. The potential involvement of the brain and its plastic mechanisms in glaucoma suggests that the diagnosis should involve the assessment of neuronal function and the treatment should consider the entire visual pathway [55][56][57] .
Evaluating the presence of neuroplasticity requires accurate experimental conditions (e.g. using matched SS). While we attempted to match the visual input between participants with glaucoma and their matched control participants in the most accurate way possible, there are limitations associated with the simulation of the VFDs. First, the simulations of the scotomas were designed based on SAP contrast sensitivity. While SAP is the gold standard to clinically access the visual field available, SAP may not accurately represent the perception of glaucoma patients (for example, it does not assess perceptual masking). Moreover, individual performance is, amongst others, affected by attention and experience 58 . Errors in the contrast sensitivity assessment will lead to inaccurate simulations and biases that, in turn, may erroneously be interpreted as signs of reorganization. Second, www.nature.com/scientificreports/ simulations primarily result in short-term plastic changes in pRFs, in particular those associated with perceptual masking 37 and spatial attention 59 . This may exacerbate the difference between glaucoma and controls and affect the assessment of long-term plastic changes associated with glaucoma. Third, the use of simulations in pRF mapping might introduce analytical artifacts (modeling biases), in particular at low SNR. Such pRF modeling biases might have the opposite effect of approximating the visual perception (during retinotopic mapping) of control participants to those of glaucoma. To test whether the simulations increased the disparity between controls and glaucoma participants, we calculated the local deviants of pRF properties between participants with glaucoma and controls without simulations (Fig. S11). Similar to the analysis based on the controls with simulations, local variations in the pRF characteristics are observed in all glaucoma participants (size, eccentricity and polar angle are shown in Figs. S12, S13 and S14, respectively). Despite these limitations, we conclude that simulations are a useful tool when assessing visual plasticity. They predict the perception of a glaucoma patient if the consequence of their disease would be a reduced retinal input only. Moreover, for various reasons, the simulations are unlikely to introduce strong biases in the pRF estimates. Firstly, the visual field reconstructions show that, although the simulations do not capture the total extent of the VF loss, they better mimic the glaucomatous VF than do controls without a simulation. Secondly, the simulations resulted in an attenuation of stimulus contrast. In a previous study, we reduced retinotopic stimulus contrast to 2% (from the conventional 50%). And while this reduced the BOLD modulation, it did not affect the pRF estimates in early visual areas 36 . Thirdly, the pRF deviants obtained with and without simulations are very similar, in particular local pRF deviations were found in all glaucoma participants. This supports that glaucoma is associated with plasticity, albeit limited in spatial extent and that little to no structural differences are induced by the SS. The participants with glaucoma were heterogeneous in various aspects, for instance in the extent of their scotoma, their disease duration, and in the asymmetry of the VFDs between both eyes. Such differences may affect the degree of neural reorganization as well as the degree of predictive masking taking place. Inside the scanner, only a relatively limited central part of the VF could be stimulated. Yet, in glaucoma, VFDs originate in the periphery of the visual field and many of the foveal scotomas that we could assess had relatively little reductions in contrast sensitivity. Nevertheless, we found marked differences in BOLD response, amongst others. Still, for this reason, future studies could consider including participants with more advanced glaucoma that would also have binocularly overlapping scotomas in their central vision. This could help to further establish the relation between the severity of the disease, the magnitude of the BOLD signal and any deviations in pRF properties. This information would also be beneficial when assessing the accuracy of the VF reconstructions. Alternatively or additionally, it could be useful to perform the studies using visual stimulation that could reach deeper into the visual periphery 60 .
The specific origin of the pRF reorganization patterns is not known nor are the factors that determine which specific adaptations take place (position shifts or size changes), this implies that some major points need to be considered. First, the techniques that we use reflect the aggregate RF properties at the population and/or subpopulation level. The pRF dynamics that we measure could result from changes in a subset of neurons 17 . Second, changes in pRF properties may also result from extra-classical RF modulations and from attentional modulation. Such interactions can be studied by applying advanced neural computational models which have the ability to capture the activity of multiple subpopulations 31,36 to take into account the extra classical RF representation 61 ; to model the effect of attention and higher order cognitive functions 59 and using stimuli that target specific neural populations 31,36 . Stimulus-driven approaches, such as pRF mapping, have an inherent disadvantage when applied to study neuroplasticity after VF loss. As we have seen in our present study, despite adding accurate control conditions, differences in the visual input of participants with glaucoma and control participants may still potentially influence results 62 . One way to circumvent this limitation is to apply cortico-cortical models which are stimulusagnostic by design. When applied to resting state data an approach such as connective field modeling may be used to evaluate if similar patterns of reorganization can be found in the absence of stimulation 63,64 .
Although this study sheds light on the neural mechanisms underlying predictive masking of natural scotomas and cortical reorganization, the entity and mechanism of cortical reorganization are still not completely clear. As considerable reorganization can be expressed subcortically, this should be studied using ultra high fields fMRI where different layers across the cortex can be measured. This will be important to quantify the level of adult brain plasticity in visual processing. Moreover, the VF reconstruction based on fMRI will improve using higher spatial resolution scans.

Conclusion
Participants with glaucoma exhibit individually unique patterns of reorganization that, nevertheless, suggest that the primary visual cortex of adults with glaucoma retain a local and limited degree of reorganization. This manifested itself in shifts in the centers of pRFs as well as in changes in their sizes. For some participants with glaucoma these changes even extend beyond their scotoma projection zone. Such changes to the neuronal configuration of the RFs may contribute to the masking of VFD which prevents patients from noticing their VFD. Moreover, although limited in spatial extent, this neuroplasticity may be critical to the successful implementation of future restorative therapies, e.g. those based on stem-cells.

Methods
Study population. Nineteen individuals with POAG and nineteen control participants with normal or corrected vision were recruited. The population demographics are shown in Table 1. The monocular data acquired for these participants was previously included in another study 32 , so the details about the inclusion and exclusion criteria as well as the details of the ophthalmic data acquisition can also be found at 32  Inclusion criteria for the participants with glaucoma were as follows: having an intraocular pressure (IOP) > 21 mmHg before treatment onset, the presence of a VFD due to glaucoma (glaucoma hemifield test outside normal limits), abnormal optical coherence tomography (OCT); peripapillary retinal nerve fiber layer thickness (pRNFL) at least one clock hour with a probability < 0.01, spherical equivalent refraction within ± 3 D.
Exclusion criteria for both groups were: having any ophthalmic disorder affecting visual acuity or VF (other than POAG in the participants with glaucoma), any neurological or psychiatric disorders, the presence of gross abnormalities or lesions on their magnetic resonance imaging (MRI) scans, or having any contraindication for MRI (e.g., having a pacemaker or being claustrophobic).
Ophthalmic data. Prior to their participation in the MRI experiments, we assessed for all participants their visual acuity, IOP, VF sensitivity (measured using HFA and frequency doubling technology [FDT]) and macula nerve fiber layer (MNFL) thickness. Visual acuity was measured using a Snellen chart with optimal correction provided for the viewing distance. IOP was measured using a Tonoref noncontact tonometer (Nidek, Hiroishi, Japan). The VFs were first screened using FDT (Carl Zeiss Meditec) using the C20-1 screening mode. The contrast sensitivity at several locations of the VF was measured using SAP specifically using a HFA (Carl Zeiss Meditec, Jena, Germany) with the 24-2 or 30-2 grid and the Swedish Interactive Threshold Algorithm (SITA) Fast. Only reliable HFA tests were included in this study. A VF test result was considered unreliable if false-positive errors exceeded 10% or fixation losses exceeded 20% and false-negative errors exceeded 10% 65 . Finally, the RNFL thickness was measured by means of OCT using a Canon OCT-HS100 scanner (Canon, Tokyo, Japan). Experimental procedure. Each participant completed two (f)MRI sessions of approximately 1.5 h each. In the first session, the anatomical scan (T1w), Diffusion Weighted Imaging (DWI), T2w, resting state functional scans and a MT localizer were acquired. In the second session, the retinotopic mapping and scotoma localizers experiments took place. These experiments were performed binocularly and monocularly as well. The resting state fMRI results are reported on in a different paper 66 . Here, we report on the results of the binocular retinotopic mapping scans.
Participants with glaucoma. The second session differed for participants with glaucoma and control participants. For the participants with glaucoma, their second (f)MRI session comprised the retinotopy and scotoma localizer experiments. The retinotopy experiment comprised nine runs in total of which six were done with binocular and three with monocular vision. For this study, only the binocular runs were analyzed. The scotoma localizer experiment comprised 2 runs in total of which one was done with binocular and one with monocular vision. This task was performed to control residual activity within the scotoma projection zone. In the monocular experiments, the most lesioned eye was stimulated and the other was occluded using an MRI compatible opaque lense. The most lesioned eye was selected based on the SAP MD (mean deviation) score; the eye with the lowest MD was selected. The monocular retinotopy results were used to assess the capability of fMRI to detect visual field defects and are reported on in a different paper 32 .

Control participants.
In their second (f)MRI session, the control participants performed the LCR, LCR SS and scotoma localizer experiments. The latter was used to define the simulated scotoma projection zone. All experiments were done with binocular vision. For both LCR and LCR SS, four runs were performed. Two scotoma localizers were acquired, one with and another without the SS superimposed on the stimulus.

Stimulus presentation and image acquisition.
Stimuli were presented on an MR compatible display screen (BOLDscreen 24 LCD; Cambridge Research Systems, Cambridge, UK). The screen was located at the head-end of the MRI scanner. Participants viewed the screen through a tilted mirror attached to the head coil. Distance from the participant's eyes to the display (measured through the mirror) was 120 cm. Screen size was Table 1. Demographics of participants with glaucoma and controls. Average and standard deviation of age, percentage of female, intraocular pressure (IOP) for the right and left eye (average over three measurements), peripapillary retinal nerve fiber layer (pRNFL) thickness and the VF mean deviation (VFMD) measured with SAP for the right and left eye. Note that participants with glaucoma were receiving treatment.

Average (n = 19) Standard Deviation (n = 19) Average (n = 19) Standard Deviation (n = 19)
Age ( www.nature.com/scientificreports/ Stimuli. All participants underwent binocular visual field mapping using luminance contrast retinotopy (LCR) mapping. Figure 7A shows an example frame of the stimulus. Additionally, the glaucoma participants observed the LCR monocularly, and the healthy participants viewed the LCR binocularly with a simulated scotoma (LCR SS) superimposed (Fig. 7B). For each control, the LCR SS was matched to that of a participant with glaucoma (see "Luminance-contrast defined retinotopy with simulated scotomas (LCR SS)"). This condition acted as a reference for the glaucoma binocular LCR, in order to disentangle possible (cortical) plasticity.
Luminance-contrast retinotopy (LCR). LCR consisted of a drifting bar aperture defined by high-contrast flickering texture 69 . The bar aperture, i.e. alternating rows of high-contrast luminance checks drifting in opposite directions, moved in eight different directions: four bar orientations (horizontal, vertical, and the two diagonal orientations) and for each orientation two opposite drift directions. The bar moved across the screen in 16 equally spaced steps, each lasting 1 TR (repetition time, time between two MRI volume acquisitions). The bar contrast, width, and spatial frequency were 100%, 1.75°, and 0.5 cycles per degree, respectively. After each pass, during which the bar moved across the entire screen during 24 s, the bar moved across half of the screen for 12 s, followed by a blank full screen stimulus at mean luminance for 12 s as well.
Luminance-contrast defined retinotopy with simulated scotomas (LCR SS). LCR SS consisted of the LCR stimulus with a simulated scotoma. The SS for a control participant was designed to mimic the contrast sensitivity of the corresponding glaucoma participant under binocular vision. The scotoma was simulated by means of local reductions in stimulus contrast. In particular, the SS consisted of an alpha transparency contrast layer defined using the HFA sensitivity values of the respective participant with glaucoma. For example, a decrease of 3 dB in www.nature.com/scientificreports/ HFA sensitivity is simulated by means of a reduction in stimulus contrast of 50%. The binocular HFA sensitivity at every position measure was calculated by taking the maximum between Left and Right eye.
Attentional task. During the retinotopic mapping scans, participants were required to perform a fixation task in which they had to press a button each time the fixation cross changed colour between black and yellow (retinotopic experiments) and between white and black (scotoma localizer). The fixation cross extended towards the edges of the screen so that it could be used as a cue for the screen's center by the participants with central scotomas. The average performance was above 75% for all conditions and for participants with glaucoma and control participants. The task performance per condition is presented in Table S1.
Magnetic resonance imaging. Data acquisition and preprocessing. Scanning was carried out on a 3 Tesla Siemens Prisma MR-scanner using an 64-channel receiving head coil. A T1-weighted scan (voxel size, 1mm 3 ; matrix size, 256 × 256 × 256) covering the whole-brain was recorded to chart each participant's cortical anatomy.
Padding was applied to strike a balance between participant comfort and a reduction of possible head motion. The retinotopic scans were collected using standard EPI sequence (TR, 1500 ms; TE, 30 ms; voxel size, 3mm 3  The T1-weighted whole-brain anatomical images were reoriented in AC-PC space. The resulting anatomical image was initially automatically segmented using Freesurfer 70 and subsequently edited manually. The cortical surface was reconstructed at the gray/white matter boundary and rendered as a smoothed 3D mesh 71 . The functional scans were analysed in the mrVista software package for MATLAB (available at https:// web. stanf ord. edu/ group/ vista/ cgi-bin/ wiki/ index. php/ MrVis ta). Head movement between and within functional scans were corrected 72 . The functional scans were averaged and coregistered to the anatomical scan 72 , and www.nature.com/scientificreports/ interpolated to a 1 mm isotropic resolution. Drift correction was performed by detrending the BOLD time series with a discrete cosine transform filter with a cutoff frequency of 0.001 Hz. In order to avoid possible saturation effects, the first 8 images were discarded.
Visual field mapping and ROI definition. The pRF analysis was performed using both conventional population receptive field (pRF) mapping 69 and micro-probing 73 . Using both models, for all the participants the functional responses to binocular LCR were analysed using a full field (FF) model (Fig. 7C). Additionally, the data acquired in LCR SS condition in the control participants in experiments 1 and 2 were analyzed using a model that included the simulated scotoma (scotoma field; SF, Fig. 7D). We observed no significant difference between the SF and FF models applied to controls SS, for this reason the results presented in this study were calculated using the FF model. The pRF models obtained for glaucoma participants were calculated uniquely using the FF model, due  www.nature.com/scientificreports/ Figure 6. Comparison of the V1 visual field reconstruction between glaucoma and control SS participants. VF reconstruction glaucoma-control participants SS pairs and the SAP tests for the left and right eye are shown. Note that to compare the binocular reconstructed VF with the monocular HFA outcomes, the highest contrast sensitivity from both eyes should be considered. The visual field reconstructions were obtained using pRFs estimated using the micro-probing technique. www.nature.com/scientificreports/ to the fact that contrary to the SS, we do not know the 'ground truth' (exact shape of the scotoma of glaucoma participants). Moreover the Micro-probing and conventional pRF mapping resulted in similar pRF estimates, the pRF analysis of eccentricity, polar angle and size (shown in the Figs. 1, 2, 3 and 4), were obtained using the conventional pRF technique, while the visual field reconstructions were estimated using micro-probing. The visual areas V1, V2 and V3 were defined on the basis of phase reversal on the inflated cortical surface, obtained with the conventional pRF model using the LCR stimulus presented binocularly.
Population receptive field analysis. As in previous work 15,69,74 data was thresholded by retaining the pRF models that explained at least 0.15 of the variance in the BOLD response and that had an eccentricity in the range of 0°-7°, for all conditions (i.e., LCR and LCR SS).
Correlation analysis. The correlations between the BOLD modulation, the pRF size and position deviation and the disease severity were calculated using a linear mixed effects model with a slope and intercept per subject as a random effect.
where y is the dependent variable, i.e. BOLD modulation, pRF size and position deviation, x is the independent variable, i.e. contrast sensitivity and MNFL, and S is the subject. To determine whether the BOLD modulation, defined as the standard deviation of the BOLD signal, of glaucoma participants was different from control participants, repeated measures two-way analysis of variance (ANOVA) with ROIs, eccentricity and condition (LCR and LCR SS) were performed.

Analysis of the variation in pRF properties between glaucoma and control participants.
To investigate changes in pRF sizes between glaucoma and control participants the cumulative distribution of pRF size is depicted across voxels. At the group level, statistical testing was performed by calculating the median pRF size across voxels per participant and comparing these median values across groups (glaucoma vs control participants) using a twosample t-test. The BOLD modulation and pRF estimates (eccentricity, polar angle and size) were binned over eccentricity, in 1° bins, i.e. the average BOLD modulation of the voxels with an estimated eccentricity between each bin, (e.g. 1°-2°), were calculated and plotted as function as the fixed (categorical) eccentricity value (e.g. 1.5°). Unless otherwise specified the error bars correspond to the confidence intervals (CI), the 95 th percentile of the full distribution of the BOLD modulation and pRF estimates aggregated across participants.
Given the heterogeneity of the VFD of the glaucoma participants, one of the challenges is how to properly perform a group level analysis. Analyzing the pRF deviation at the level of quarter fields and correlating these deviations with SAP and OCT metrics, allows us to understand if at the group level, the deviations in pRF positions are related to disease severity. In order to understand how pRF properties are associated with the severity of glaucoma, we performed two different analysis: (1) we correlated the averaged pRF size and ED of the average position calculated per quarter field between participants with glaucoma and respective control participant with the contrast sensitivity and MNFL thickness; and (2) we assessed the significance of the deviation between pairs in relation to an expected baseline deviation, which summarizes the variation in pRF properties amongst control participants. In order to take this into account, a group level analysis was done by determining the rank of the deviation of the participant pair (either the glaucoma-control SS pair or the control SS-control NS pair) relative to the baseline deviation (which consisted of all the 19 deviations between the matching control participant and all other control participants, Fig. S10). This approach and its results are presented in detail in the Sects. 7 and 8 of the Supplementary Information.
Visual field reconstruction analysis. Both MP and pRF mapping techniques allow an accurate reconstruction of the VF by back-projecting the pRF properties of all voxels within a visual area onto the VF 32,33 . The VF reconstructed maps reflect the VF sampling density. Since the presence of VFD reduces the sampling of a particular region of the VF, fMRI-based VF reconstruction techniques are suitable to detect VFDs and indirectly reflect VF sensitivity.
In this study we compared the visual field reconstruction between participants with Glaucoma and Controls SS. The visual field reconstruction was obtained using the MP technique as described in 32 . In addition, to directly compare the VF reconstructions between groups, VF maps were converted from normalized scale to a dB scale by taking the 10 × log10 of the sampling density values, resulting in VF sensitivity values. The final VF maps correspond to the deviation of VF sensitivities between glaucoma and controls SS from controls NS according to the 90, 95, 98, 99 and 99.5% CI boundaries. This approach was previously applied by 33 .
Statistical analysis. All statistical analyses were performed using R (version 0.99.903; R Foundation for Statistical Computing, Vienna, Austria) and MATLAB (version 2016b; Mathworks, Natick, MA, USA). After correction for multiple comparisons, a p-value of 0.05 or less was considered statistically significant.

Data availability
The raw data supporting the conclusions (ophthalmic evaluations and MRI data) of this article are publicly available at Dataverse NL under the project name 'Local neuroplasticity in adult glaucomatous visual cortex' , https:// doi. org/ 10. 34894/ QMADIQ. The code for the Micro Probing based mapping of RFs and visual field reconstruction is available via the following links: https:// github. com/ Joana-Carva lho/ Micro-Probi ng; https:// github. com/ Joana-Carva lho/ FMRI-based-Visual-field-Recon struc tion.