Analysis of correlations between local geographic atrophy growth rates and local OCT angiography-measured choriocapillaris flow deficits

: The purpose of this study is to quantitatively assess correlations between local geographic atrophy (GA) growth rates and local optical coherence tomography angiography (OCTA)-measured choriocapillaris (CC) flow deficits. Thirty-eight eyes from 27 patients with GA secondary to age-related macular degeneration (AMD) were imaged with a commercial 1050 nm swept-source OCTA instrument at 3 visits, each separated by ∼ 6 months. Pearson correlations were computed between local GA growth rates, estimated using a biophysical GA growth model, and local OCTA CC flow deficit percentages measured along the GA margins of the baseline visits. The p-values associated with the null hypothesis of no Pearson correlation were estimated using a Monte Carlo permutation scheme that incorporates the effects of spatial autocorrelation. The null hypothesis (Pearson’s ρ = 0) was rejected at a Benjamini-Hochberg false discovery rate of 0.2 in 15 of the 114 visit pairs, 11 of which exhibited positive correlations; even amongst these 11 visit pairs, correlations were modest ( r in [0.30, 0.53]). The presented framework appears well suited to evaluating other potential imaging biomarkers of local GA growth rates.

Recently, there has been interest in associating OCTA-based metrics of CC impairment with GA growth rates [11][12][13][14][15][16]. Thus far, studies have focused primarily on assessing correlations between CC impairments and global measures of GA growth (i.e., measures that do not resolve how different segments of a lesion margin expand) [11][12][13][14]16]. These studies have reported positive correlations between CC flow deficits (FDs) and global GA growth rates, but disagree about whether CC FDs in regions adjacent to GA margins exhibit different correlations than do CC FDs in regions farther from the margins. To our knowledge, there have thus far been two studies of local (i.e., spatially resolved along a lesion margin) correlations between CC impairments and GA growth rates [11,15]. In a recent study by Sacconi et al., GA lesions were reported to preferentially expand into regions of CC impairment [15]. In contrast, in a pilot study by our group, we found no obvious correlation between local CC impairments and local GA growth rates [11]. However, the primary purpose of that study was to demonstrate the applicability of our analysis framework, and its small enrollment precluded a rigorous statistical analysis.
In this study we apply our previously developed framework to assess correlations between local CC FDs and local GA growth rates in a larger cohort of GA eyes longitudinally imaged at 3 visits with a commercial 1050 nm swept-source OCTA (SS-OCTA) instrument. Additionally, we apply spatial statistical methods to quantify the statistical significance of the observed correlations. A central aim of this study is to demonstrate the utility of our framework for rigorous statistical assessment of correlations between a local imaging biomarker (here, CC impairment) and local GA growth rates. From a pathogenesis viewpoint, such a framework is interesting because, noting that GA growth is highly anisotropic, if an imaging biomarker were a surrogate for a dominant, causative driver of local GA growth, we would expect there to be a local correlation. From the viewpoint of GA therapeutic trials, such a framework is important because if a locally correlated imaging biomarker were identified, it could be used to predict GA growth patterns, providing a more complete patient stratification compared to global imaging biomarkers.

Patient enrollment and imaging of GA and CC FDs
Patients with GA secondary to non-exudative AMD were enrolled in a prospective OCTA imaging study that was approved by the institutional review board of the University of Miami Miller School of Medicine. Enrollment was from June 2016 through November 2019. Informed consent was obtained from each subject. This study was performed in accordance with the tenets of the Declaration of Helsinki and complied with the Health Insurance Portability and Accountability Act of 1996. To be included in the study, eyes were required to: (1) have underwent SS-OCTA imaging at a baseline visit (Visit 1), a 6 month follow-up visit (Visit 2), and a 1 year follow-up visit (Visit 3); (2) have a total GA area ≥2.54 mm 2 (one disc area), and, for multifocal lesions, have at least one GA focus with an area ≥1.25 mm 2 ; and, (3) have all GA foci fully contained within the 6 mm × 6 mm field-of-view at all visits. Additionally, eyes were excluded if the GA was continuous with parapapillary atrophy, if there was macular atrophy due to a diagnosis other than non-exudative AMD, if there was any history of exudative macular neovascularization (MNV), or if treatment naïve, non-exudative MNV was identified by SS-OCTA imaging [17].
OCT and OCTA imaging was performed using a commercial SS-OCTA instrument (PLEX Elite 9000; Carl Zeiss Meditec, Dublin, CA), as previously described [12]. This SS-OCTA instrument has a swept laser light source with a central wavelength of 1050 nm, a bandwidth of 100 nm, and full-width-at-half-maximum axial and transverse resolutions of approximately 5 µm and 20 µm in tissue, respectively. The instrument operates at a speed of 100,000 A-scans per second. The 6 mm × 6 mm scan pattern, which was used for both OCT and OCTA imaging in this study, consisted of 500 A-scans per B-scan along the 6 mm transverse dimension, with each B-scan repeated 2 times to generate the OCTA signal, and 500 B-scan positions along the 6 mm vertical dimension, resulting in an isotropic 12 µm A-scan spacing. The interscan time (i.e., the time between repeated B-scans) was 6.1 ms. Each A-scan had a depth of 3 mm in tissue, and was comprised of 1536 pixels. Images with a signal strength of <7 were excluded from analysis, as were images with pronounced motion artifacts.
En face structural OCT images were used to identify GA on the basis of OCT hyper-transmission on a slab with segmentation boundaries spanning 64 µm to 400 µm below Bruch's membrane [12,18], and the presence of GA was confirmed on corresponding B-scans per the Consensus on Atrophy Meetings classification [3]. After an eye was enrolled into the study based on the minimum size requirements for unifocal and multifocal lesions, all additional areas of hypertransmission in the scan area having a greatest linear dimension of at least 250 µm were included in the overall lesion measurements. The boundaries of GA were manually outlined by two independent graders (YS and LW) using image analysis software (Adobe Photoshop CC; Adobe Systems, San Jose, CA) and a consensus outline was reached by both graders. If a consensus could not be reached, then a senior grader (PJR) served as the adjudicator. Inconsistencies in GA tracings between visits were managed by imposing a non-shrinkage condition, whereby GA tracings from earlier visits were clipped to be inside of, or coincident with, GA tracings from later visits. Regions of tracing inconsistency-that is, image regions between the clipped and un-clipped lesion contours-were excluded from subsequent CC OCTA analysis.
En face OCTA CC FD images were formed for Visit 1 and Visit 2 acquisitions as described in our previous studies [12,19,20]. Briefly, CC blood flow was detected using the complex optical microangiography (OMAG c ) algorithm [6,[21][22][23]. En face CC OCTA images were computed by projecting the OCTA volume through a 16 µm thick slab having an inner boundary positioned 4 µm posterior to Bruch's membrane, which was segmented using a semi-automated software; the accuracy of this segmentation was visually confirmed, and segmentations were manually adjusted where required. The CC slab thickness was selected considering the physiological thickness of the CC and the axial resolution of the instrument; the CC slab position was selected considering the anatomical position of the CC, which directly abuts Bruch's membrane [24]. In particular, a 4 µm offset from Bruch's membrane was selected so as to reduce background OCTA signal (i.e., noise) generated by Bruch's membrane [24]. After removing retinal vasculature projection artifacts, OCTA signal compensation was performed using normalized and inverted (logarithmic-scale) en face (structural) OCT images from the corresponding slab. Choriocapillaris FD images were computed via global thresholding of the en face CC OCTA images, as described elsewhere [25][26][27][28]. After thresholding, any CC FDs with a diameter smaller <24 µm, which correspond to physiologically normal CC FDs and speckle noise, were removed [29].

Visit registration
The Visit 2 and Visit 3 CC FD images and lesion tracings were mapped into the Visit 1 coordinate frame using corresponding landmarks-retinal bifurcations-manually selected from en face OCTA projections of the retinal vasculature, which were formed using an automatic segmentation. In particular, the selected landmarks were used to estimate second-order polynomial transformations [30] that aligned the Visit 2 and Visit 3 coordinate frames with the Visit 1 coordinate frame. For CC FD images, resampling with linear interpolation was performed following alignment.

Quantification of local GA growth rates
Following our previous study [11], we estimated local GA growth rates via a GA growth modeling approach. In this approach we consider GA growth as an interface propagation problem, wherein the GA margin of the baseline visit is enlarged, according to some specified dynamics, until it is coincident with the GA margin of the follow-up visit. Specifically, letting Ω ⊂ R 2 be the image domain, we denote the observed GA region and GA margin at time t by G(t) ⊂ Ω and ∂G(t) ⊂ Ω, respectively, and consider GA growth for times t ∈ [t b , t f ], where t b is the time of the baseline visit, and t f is the time of the follow-up visit. Furthermore, we denote the modeled (i.e., estimated) GA region and GA margin at time t byĜ(t) ⊂ Ω and ∂Ĝ(t) ⊂ Ω, respectively. Because GA lesions typically undergo merging, either between different segments of the same lesion focus, or between segments of different lesion foci, we make use of the level set method [31], which representsĜ(t) and ∂Ĝ(t) via level sets of the function ϕ(x, t) : Ω × R → R. In particular, ϕ is constructed such thatĜ(t) = {x ∈ Ω : ϕ(x, t) ≤ 0} and ∂Ĝ(t) = {x ∈ Ω : ϕ(x, t) = 0}. With this formulation, our boundary conditions can be expressed as: or, equivalently: where ϕ b (x) and ϕ f (x) are the signed distance functions of ∂G(t b ) and ∂G(t f ), respectively, with the convention that regions within the GA lesion are assigned negative distances. In general, we can model GA growth by evolving ϕ according to: where ∂ t is the partial derivative operator with respect to time, V determines the rate of lesion expansion, and ∇ is the gradient operator. For our current study, we adopt the simplified GA growth dynamics specified by: where α and β are positive constants (model parameters) and κ(ϕ) = ∇ · (∇ϕ/||∇ϕ||) is the curvature of ∂Ĝ(t). The parameter α causes a constant outward lesion growth, and the parameter β causes concave lesion segments to grow faster than convex lesion segments; geometric and physiological rationales for β are presented in Fig. 1. To ensure that the boundary conditions of Eq. (2) are satisfied, when solving Eq. (3) we enforce the constraint that: After solving Eq. (3), for each x b ∈ ∂G(t b ) we construct a growth trajectory, γ(x b , t), which describes the spatial path a hypothetical particle positioned at the space-time point (x b , t b ) traces to reach a point in {(x, t f ) : x ∈ ∂G(t f )}. Note that these growth trajectories are under-determined because a sequence of level sets does not uniquely determine a point correspondence [32]. To address this ambiguity, we adopt the assumption of zero tangential velocity. Rationale for this assumption is provided in Fig. 1(C). In particular, note that, for circular, isotropically growing lesions, this assumption yields growth trajectories whose lengths agree with the commonly used square-root-of-area growth rate metric, scaled by a factor of 1/ √ π [18,33] ( Fig. 1(C1)). With this assumption, growth trajectories γ(x b , t) can be constructed by integrating [32]: on the time interval t ∈ [t b , t f ] subject to the initial condition that γ(x b , t b ) = x b . As in Moult et al. [11], merging trajectories are excluded from the analysis, with the rationale that tracking becomes highly ambiguous after merging. Example growth trajectories are shown in Fig. 2(B).

GA growth model parameters and numerical implementation
The same model parameter values of α = 1 and β = 0.75 that have been used in our previous growth model studies [11,34] were used in the current study. These model parameters were Note that the trajectory is orthogonal to the tangent line at p (indicated by the dashed blue line segment), meaning that its component tangential to the lesion boundary, γ | | , is zero. This produces growth trajectories that are coincident with the circles' radii vectors, and whose lengths agree with the the commonly used square-root-of-area growth rate metric, scaled by a factor of 1/ √ π [18,33]. Non-normal growth results in growth trajectories such as the one shown in purple, which have a non-zero tangential component. (C.2) Microscopic interpretation: under reasonable conditions, progressively 'zooming' in on a point p on any GA margin yields an approximately straight line. Ignoring the curvature term in Eq. (4) (i.e., β = 0), and letting the time-interval of growth become small, the normal growth assumption results the growth trajectory shown in red, whose component tangential to the lesion boundary, δγ | | , is again zero. Scatterplot showing local GA growth rates as a function of local CC FD percentages, with each marker corresponding to a unique position on the Visit 1 GA margin. The points within the dashed ellipse correspond to a contiguous, high-growth segment of the GA margin (see corresponding dashed ellipses in panels B and C). The sample Pearson correlation coefficient, r, between the local GA growth rate and local CC FD percentage can be computed from the scatterplot measurements; however, because the samples are autocorrelated, additional analyses must be performed to compute the associated p-value (see Fig. 3).
subjectively selected for yielding qualitatively plausible growth trajectories. In general, as β increases in value relative to α, the growth trajectories become smoother (see Fig. 1(B)); however, the resulting trajectories are, for most lesion segments, only mildly perturbed by moderate changes in β.
All numerical analysis was performed in MATLAB 2019b (MathWorks, Inc., Natick, MA). The grid spacing ∆x used for the analysis matched the A-scan spacing (i.e., ∆x = 12 µm). Eq. (3) was numerically solved using the toolbox developed by Mitchell et al. [35]. In particular, the term α||∇ϕ(x, t)|| was computed using a second order essentially nonoscillatory scheme, and the term −βκ(ϕ)||∇ϕ(x, t)|| was computed using a second order central difference scheme [36]. These terms were numerically integrated using a first order forward Euler integrator, with a fixed timestep of ∆t = 0.1 (arbitrary units). The level set was reinitialized after every timestep by iteratively solving the reinitialization equation [37]. The GA growth model was evolved until steady state (i.e., when the baseline GA margin is coincident with the follow-up GA margin), and Heun's method was used to integrate Eq. (6).

Quantification of local CC impairment
Following Moult et al. [11], local CC impairments were assessed by computing the local CC FD percentages along the GA margin of the baseline visit (i.e., Visit 1 for the Visit 1-Visit 2 pair; Visit 2 for the Visit 2-Visit 3 pair; and Visit 1 for the Visit 1-Visit 3 pair). Specifically, for each measurement position p on the GA margin, the CC FD percentage was computed within a 500 µm diameter disk template centered on p (Figs. 2(E), 2(F)). Lesion interiors were excluded from the computation; areas corresponding to tracing inconsistencies between visits, described previously, were treated as missing data and were also excluded from the analysis. Measurement positions were located at 6 µm intervals (in arclength) along the baseline GA margin.

Statistical assessment of local correlations via variogram matching
To assess the strengths of correlations between local CC impairments and local GA growth rates we used statistical hypothesis testing ideas: for each visit pair (i.e., Visit 1-Visit 2, Visit 2-Visit 3, and Visit 1-Visit 3), we tested the null hypothesis that local GA growth rates and local CC FD percentages are Pearson uncorrelated (i.e., Pearson's ρ = 0). As described in previous sections, our analysis yields a set of paired local CC FD percentage measurements and local GA growth rate measurements at 6 µm intervals along the baseline GA margin. Standard procedures for testing ρ = 0 under independence assumptions are not applicable in our study because the measurements are highly autocorrelated. To this end, following Viladomat et al. [38] and our prior study [34], we estimate a suitable null distribution that considers spatial autocorrelation (Fig. 3). As illustrated in Fig. 3, the basic strategy of this approach is to randomly permute the spatial positions of the local CC FD percentage measurements, and then to spatially smooth the randomly permuted values in such a way that their spatial autocorrelation matches that of the original, unpermuted data. This smoothing step involves matching the variograms of the permuted CC FD percentage measurements to the variogram of the original, unpermuted CC FD percentage measurements (see Section 2.7 for details). To manage the computation cost of estimating the null distribution, for lesions with more than N s = 1000 margin positions, a random sample of size N s was selected.

Variogram matching parameters and numerical implementation
Our p-value computation for correlations between local GA growth rates and local CC flow deficits follows nearly identically the steps described in Viladomat et al. [38]. We do, however, diverge in our selection of the optimal smoothing bandwidth, δ (See Appendix B of Viladomat et al. for details of this parameter). Specifically, with the aim of decreasing computation cost, rather than computing δ for each random permutation (as in Viladomat et al.), we compute a single The CC FD percentage measurements are randomly permuted, creating a sample from a random field that is independent of the random field from which the local GA growth rate measurements came. The spatial structure of these randomly permuted measurements is then matched to the spatial structure of the unpermuted measurements by variogram matching. (D) The Pearson correlation coefficient, r, between the permuted and variogram-matched local CC FD percentage measurements and the unpermuted local GA growth rate measurements is computed. This process is repeated N MC = 10,000 times to construct a null distribution (light gray histogram) against which the observed correlation may be validly compared. The null distribution under the assumption of independent samples (dark gray histogram) has a smaller variance compared to that of the derived null distribution that accounts for the spatial autocorrelation of the samples. In this particular visit pair, the Pearson correlation for the observed (unpermuted) data is r = -0.30 (see Fig. 2(G)), which corresponds to a vanishingly small p-value under the assumption of independence, but a much larger p-value when spatial autocorrelation is incorporated in the calculation.
value of δ using a multiresolution search on a small, initial set of random permutations, and apply this selected δ value for all N MC random permutations used to estimate the null distribution (Fig. 4). In particular, we first perform a 'coarse' parameter search over a parameter set ∆ c , computing the mean squared error between the target (observed) variogram Λ T of the local CC FD percentage, and M variograms corresponding to different random permutations of the local CC FD percentage: where ∆ c ∈ {ϵ c , 2ϵ c , . . . , 1 − ϵ c }, ϵ c is a positive constant, which determines the coarseness of the search, π m is the m-th random permutation of the measurement values, Λ π m its associated variogram, and L = {0, 12, . . . , 12 × ⌊ℓ 25 th /12⌋} [µm] is the set of spatial lags, with ℓ 25 th being the lag corresponding to the 25 th percentile of the distribution of measurement pair distances.
In the final step, we perform a 'fine' parameter search over ∆ f , using a new set of K random permutations of the local CC FD percentages: f <ϵ c again determines the coarseness of the search, and we assume for convenience that ϵ c is an integer multiple of ϵ f . The smoothing parameter δ * f is then applied for all N MC random permutations used to estimate For each eye, the optimal smoothing parameter, δ, was selected by a multi-resolution search wherein we minimized the mean squared error between the target variogram of the observed local choriocapillaris (CC) flow deficit (FD) percentages and the variograms corresponding to a set of simulated local CC FD percentages. Note that the latter correspond to random permutations of the former (see Fig. 3). After an initial coarse search with parameter spacing ϵ c = 0.1, the search was repeated over a reduced search space with a fine parameter spacing ϵ f = 0.01. The minimum argument of this second search, whose value is indicated by the red asterisk, was taken as the smoothing parameter to be used for variogram matching for the N MC = 10,000 Monte Carlo runs used to estimate the null distribution (see  Fig. 4(B).

Additional statistical analysis
After p-values were computed for each visit pair, the Benjamini-Hochberg procedure [39] was applied to obtain a candidate set of visit pairs for which we rejected the null hypothesis at a false discovery rate (FDR) of 0.2 (MATLAB function mafdr with BHFD flag set to true). For a qualitative comparison of the sample correlation coefficient distributions across different visit pairs, violin plots were generated using the custom MATLAB function by Hoffmann [40]. This function uses the built-in MATLAB function ksdensity, which was also used to estimate the optimal smoothing bandwidth for the violin plots.

Results
Data from 38 eyes from 27 patients were included in the study. A subset of this data has been used in other studies from our group, which assessed correlations between global and zonal CC FDs and global GA growth rates: Shi et al. [16] used the Visit 1-Visit 3 data from the same 38 eyes, and Thulliez et al. [12] used the Visit 1-Visit 3 data from a 22-eye subset. The mean ± standard deviation inter-visit time was 0.51 ± 0.06 years for Visit 1-Visit 2, 0.50 ± 0.05 years for Visit 2-Visit 3, and 1.01 ± 0.06 years for Visit 1-Visit 3. The mean ± standard deviation area removed by lesion clipping (as a result of tracing inconsistencies), expressed as a percentage of the baseline visit, was 4.3 ± 3.7% for Visit 1 and 2.7 ± 1.8% for Visit 2 (note that Visit 3 lesions were not clipped; see Methods). For ease of reference, eyes were labelled Eye 1 through Eye 38, arranged by increasing Visit 1-Visit 3 square-root-of-area growth rate [18,33]. Example analyses for Eye 26 (the strongest correlation) and Eye 21 (a more typical correlation) are provided in Fig. 5 and Fig. 6, respectively. The sample Pearson correlation coefficients and corresponding p-values for all 114 visit pairs are graphically summarized in Fig. 7. Applying the Benjamini-Hochberg procedure to control the FDR at 0.2, the null hypothesis was rejected in 15 of the 114 visit pairs; for ease of reference, we denote this group of 15 visit pairs by S 1 and the denote the other 99 visit pairs, which were not rejected at an FDR of 0.2, by S 0 . Within group S 1 , 4 visit pairs had negative correlations (we refer to these 4 visit pairs by S − 1 ) and 11 visit pairs had positive correlations and (we refer to these 11 visit pairs by S + 1 ). Even amongst these 11 visit pairs, correlations were modest (r in [0.30, 0.53]).

Discussion
This study used commercial 1050 nm SS-OCTA imaging in combination with GA growth modelling and spatial statistical methods to investigate correlations between local OCTAmeasured CC FDs and local GA growth rates in 38 eyes imaged at 3 visits, with each visit separated by ∼6 months. The null hypothesis (Pearson's ρ = 0) was rejected at a Benjamini-Hochberg FDR of 0.2 in 15 of the 114 visit pairs. At an FDR of 0.2, we would expect that, on average, 3 of these 15 visit pairs are false positives. While our findings do not provide support for a strong and robust correlation between local, OCTA-measured CC FDs and local GA growth rates, we emphasize that they in no way preclude the existence of a correlation between local physiological CC flow impairments and local GA growth rates; potential reasons for why there may be discrepancies between measured and physiological CC flow impairments are explored towards the end of the Discussion.
A question naturally arising from our results is whether the eyes and/or patients of group S + 1 have differing characteristics from those of group S 0 and S − 1 (see Results for the definitions of these subgroups). To explore this question, we performed a post hoc analysis of several covariates (Fig. 8). Qualitatively, for all of the examined covariates, there was substantial overlap in the covariate values from the visit pairs in S 0 and S − 1 and those from the visit pairs in S + 1 . For each of the covariates, we estimated p-values via a constrained Monte Carlo permutation test under the null hypothesis that the mean of the covariate values of S 0 and S − 1 equaled that of S + 1 (see Supplement 1 for testing details; p-values are provided in Fig. 8). This covariate analysis should be interpreted with care because S + 1 is comprised of a small subset of visit pairs, which come from an even smaller subset of eyes and subjects. Indeed, 4 of the visit pairs of S + 1 correspond to a single subject (Eye 24 [OS] and Eye 26 [OD]). As a topic of future investigation, it may be interesting to expand our search for explanatory covariates, using additional imaging data, demographic information, and/or genotype analysis. However, the clinical relevance of such investigations is likely tempered by the apparent infrequency of such subjects, and by the modest correlations. For a potential GA growth biomarker such as CC FD percentage, it seems reasonable to consider at what time scale that biomarker might be most correlated with local GA growth. In this study, there were 6 visit pairs in S + 1 corresponding to 6-month intervals (i.e., Visit 1-Visit 2 or Visit 2-Visit 3) and 5 visit pairs in S + 1 corresponding to 12-month intervals (i.e., Visit 1-Visit 3). Taken together with qualitative inspection of the scatterplots, box plots, and violin plots, which also showed similarities amongst the three visit pair types, we do not believe that this study provides evidence that imaging at 6-month versus 1-year intervals yields meaningfully different correlations between local CC impairments and local GA growth rates. However, because the temporal dynamics of CC flow-and CC flow impairment-are not fully understood, we thought it would be interesting to investigate whether some combination of CC FDs measured at different timepoints might yield stronger correlations, due to, for example, regions of dynamic CC FDs. To this end, we performed a post hoc analysis wherein we assessed the correlation between local CC FD percentages extracted from a composite CC FD image-computed by averaging the Visit 1 and Visit 2 CC FD images (Fig. 9)-with the Visit 2-Visit 3 local GA growth rates. The results of this analysis do not suggest that combining CC FD images from multiple visits provide any obvious increase in the strength of the local correlations.
In our prior pilot study [11], we used a 400 kHz A-scan rate prototype SS-OCT instrument to examine correlations between local CC impairments and local GA growth rates in 7 eyes from 5 patients. With a similar growth modeling approach as used in the current study, we qualitatively found no obvious evidence of strong local correlations. However, the small number of eyes, variable follow-up times, and lack of statistical analysis precluded drawing general conclusions. To the best of our knowledge, the study by Sacconi et al. [15], is the only other investigation of correlations between local CC impairments and local GA growth rates. In that study, which considered 30 eyes of 20 patients with a 1-year follow-up, Sacconi et al. partition the region <500 µm from the baseline GA margin into two pixel sets: one that developed atrophy at the follow-up visit (which they termed the "expansion area") and another that did not develop atrophy at the follow-up visit (which they termed the "area surrounding GA margin minus expansion area"). With this methodology, Sacconi et al. reported that the CC perfusion density in the "expansion area" was less than the CC perfusion density in the "area surrounding GA margin minus expansion area." However, comparing the CC perfusion densities in these two regions may be misleading because there is no guarantee that the distribution of the pixel-to-margin distances is the same between these two partitions. Indeed, since GA tends to expand along its margins, the pixels in the "expansion area" partition will tend to be closer to the baseline GA margin than the pixels in the "area surrounding GA margin minus expansion area" partition. Because it has been reported, both by histopathology and in vivo imaging [11,[41][42][43][44][45], that the CC is more impaired in the region immediately surrounding the lesion margin, the unequal spatial distributions of the two partitions potentially introduces a bias towards the "expansion area" pixels having more impaired CC. To illustrate this potential bias, consider the example of a unifocal circular lesion having a 500 µm radius at the baseline visit and a 750 µm radius at the follow-up visit. By the methodology of Sacconi et al., the "expansion area" would be the annulus bounded by the 500 µm and 750 µm radii circles, and the "area surrounding GA margin minus expansion area" is the annulus bounded by the 750 µm and 1 mm radii circles. Thus, all the pixels in the "area surrounding GA margin minus expansion area" are farther from the margin than are the pixels in the "expansion area," and we would expect that, due only to the analysis methodology, the CC impairment in the "expansion area" partition would be greater than that in the "area surrounding GA margin minus expansion area".This may, in part, explain the difference between the findings of our study and those of Sacconi et al. Another potential contributor to the difference between our findings is the differing CC slab positions. In particular, Sacconi et al. appear to have used an offset of the instrument's default RPE-centerline segmentation, which may not always yield a slab that is well positioned for CC FD quantification. In any case, further studies and additional analyses are needed to better understand the relationships between local OCTA-measured CC FDs and local GA growth rates.
In addition to these two prior studies of local correlations, there have been several OCTA studies investigating associations between zonal and global CC FD percentages and global GA growth rates [12][13][14]16]. Nassisi et al. quantified CC FD percentages in two concentric, 500 µm wide lesion-centered rings in a 2-visit study of 33 GA eyes [13]. They reported statistically significant associations between CC FD percentages in the first ring and global GA growth rates, as well as between the difference of the CC FD percentages in the first and second rings and global GA growth rates; however, they did not observe a statistically significant association between CC FD percentages in the second ring and global GA growth rates, nor between global CC FD percentages and global GA growth rates. In another study by the same group, Alagorie et al. quantified CC FD percentages in concentric, 100 µm wide lesion-centered rings in a 2-visit study of 30 GA eyes [14]. Correlating these CC FD percentages with global GA growth rates, they found statistically significant Spearman correlations in rings situated within 500 µm of the GA margin, but not for rings farther than 500 µm. A limitation of both of these studies is their use of 840 nm SD-OCT instruments, which, compared to 1050 nm SS-OCT instruments, have reduced penetration beneath the RPE [46][47][48]. Contrasting these studies, in a 2-visit SS-OCTA investigation of 22 GA eyes-which, as mentioned in Results, were also included in this study-Thulliez et al. reported that global CC FD percentages exhibited stronger correlations with global GA growth rates than did CC FD percentages in the 300 µm ring immediately surrounding the lesion [12]. In a follow-up study by the same group of 38 GA eyes-which were also included in this study-Shi et al. again found statistically significant correlations between global CC FDs and global GA growth rates. However, with the expanded cohort, they found no statistically significant difference in correlations between CC FDs and GA growth rates when CC FDs were considered over the entire field-of-view (i.e., global; excluding the GA region), the 0 µm − 300 µm from-the-margin-ring, the 300 µm − 600 µm from-the-margin-ring, the 0 µm − 600 µm from-the-margin-ring, or the zone comprised of the subset of the field-of-view >600 µm from the GA margin. While not directly comparable, our results seem to be most congruent with those of Thulliez et al. and Shi et al: if zonal CC impairments surrounding the lesion margin were more predictive of global GA growth rates than global CC impairments, it seems likely-although not necessary-that local correlations would also be present.
In this study, as in our prior study [34], we used the approach of Viladomat et al. to assess the statistical significance of local correlations. To better motivate our choice of this technique, it is helpful to contrast it to the standard t-test approach, which involves: (1) computing Pearson's sample correlation r between the local CC FD percentages and the corresponding local GA growth rates; (2) computing the test statistic t = r √︁ where N is the number of measurement pairs; and, (3) computing the p-value, assuming that under the null hypothesis (i.e., ρ = 0) the test statistic has a Student's t-distribution with N − 2 degrees of freedom. This approach is not applicable to our analysis because it assumes that the measurements are independent-an assumption violated by the substantial spatial autocorrelations present amongst both the local CC FD percentage measurements and the local GA growth rate measurements. That is, measurements of CC FD percentage at adjacent margin pixels are likely to have more similar values than measurements at margin pixels separated by larger distances-and similarly for the local GA growth rate measurements. This spatial autocorrelation arises both from the underlying autocorrelations in the biophysical CC impairment and GA growth processes and, in the case of CC impairment, from the finite OCT spot size and the smoothing effect of the 500 µm diameter measurement neighborhood. Statistically, the autocorrelations effectively reduce the degrees-of-freedom of the null distribution, causing the t-test approach to suffer increased Type I errors (i.e., erroneously low p-values) [38,49]. In contrast, the approach by Viladomat et al. constructs the null distribution using random fields that are independent of the field from which the measurements came and that also have matching autocorrelation structures.
It is also worth briefly mentioning our use of the Benjamini-Hochberg procedure to manage our testing of multiple hypotheses. In particular, it is reasonable to assume that our computed p-values are positively correlated at both the visit-level (i.e., p-values from visits from the same eye) and eye-level (i.e., p-values from two eyes from the same patient). However, the Benjamini-Hochberg procedure has the so-called positive regression dependency on each one from a subset (PRDS) property [50], which implies that the FDR is still controlled, even in such cases of positive dependencies.
There are several limitations to our study. First, the number of analyzed eyes is relatively modest, and therefore it is unclear how our findings generalize. This limitation is somewhat mitigated by our having OCTA data from 3 consecutive visits at regular 6-month intervals, which allows us to examine correlations at 6-month (Visit 1-Visit 2 and Visit 2-Visit 3) and 1-year (Visit 1-Visit 3) intervals. Another limitation of our dataset is that our 6-month and 1-year follow-up intervals are sufficiently long so as to make it plausible that, in some cases, new lesion foci may have both developed and merged with existing foci during the follow-up interval, which would confound the interpretation of the corresponding growth rate measurements. However, the similarity in correlations between the 6-month and 1-year imaging data, noted previously, somewhat rebuts this as a factor that substantively impacted our analysis. Our study also has methodological limitations, some specific to our study, which we shall address first and most fully, and some that are common to other OCTA CC studies, which we shall more briefly summarize. A potential limitation of our statistical analysis is that we measured correlation using Pearson's r, which is commonly used for capturing linear relationships and can be skewed by outlier measurements. Because there is no reason to expect that the relationship between local CC FD percentages and local GA growth rates is linear, and because our data clearly contain outliers (e.g., the high growth segments of Fig. 2(G)), the usage of Pearson's r may not be ideal. Our decision to use Pearson's r arose primarily from technical challenges in applying autocorrelation adjustments to rank correlations (e.g., Spearman's correlation)-in particular, we were unable to achieve well matched rank variograms. Ultimately, considering the need to control for autocorrelations and the difficultly in matching rank variograms, we opted to use Pearson's r. Despite its potential limitations, from qualitative assessment of the scatterplots of local GA growth rates versus local CC FD percentages, we do not have reason to believe that our use of Pearson's correlation-rather than, for example, Spearman's correlation-was a factor that substantively impacted our conclusions. Nevertheless, the development of techniques combining autocorrelation adjustments with other notions of correlation is an interesting future direction. Approaches such as removal or separate treatment of outliers may also be fruitful.
There are also potential limitations that arise from the variation of CC FD percentages with distance from the fovea. Specifically, comparing 3 non-overlapping, concentric rings of interest, Zheng et al. found that CC FD percentages were highest in the 1 mm fovea-centered circle, lower in the 1-3 mm annulus, and lowest in 3-5 mm annulus [20]. Moreover, these changes became more pronounced in older patients, particularly those in the 50-89-year-old range, which comprise the majority of the patients in our study. There are at least two ways that this normal variation might impact our analysis. First, the procedure of Viladomat et al. utilizes variograms, which are estimated under the assumption that the random field from which the measurements arise has a stationary (i.e., constant) mean; for the reasons above, this assumption does not hold, at least not in normal controls. Second, and related to the first, a given CC FD percentage measurement taken at a margin segment located near the fovea may not carry the same interpretation as a similarly valued CC FD percentage measurement taken at a margin segment located farther from the fovea. For example, consider a CC FD percentage value, X, measured at a margin segment located a distance R 1 from the fovea, and suppose that X is lower (i.e., less impaired) than the average CC FD percentage value,X R 1 , of all measurements from margin segments of distance R 1 from the fovea. It is plausible that when X is compared to CC FD percentage values measured at segments farther from the fovea, say at a distance R 2 >R 1 , having an average valueX R 2 , due only to physiologically normal spatial trendsX R 1 − X>X R 2 − X. Thus, by considering all margin segments together, the ability to detect X as different from other measurements along the margin can be lost, reduced, or even inverted, the latter in this case meaning that X appears to be relatively more impaired. A possible strategy that would remedy these prior two limitations is to detrend the data before analysis, effectively removing the 'normal' CC spatial variation. However, such a scheme may be challenging in eyes in which CC alterations are dominated by pathology, and in which the dependencies found in studies on normal control eyes may not directly apply. An alternative strategy would be to only correlate measurements coming from margin segments at similar distances from the fovea.
A further limitation of our analysis of CC impairment arises from the inherent ambiguity in choosing the most appropriate neighborhood for computing local CC FD percentages. In this study we opted to use a 500 µm diameter circular neighborhood, centered on the GA margin. While the neighborhood shape may be, perhaps, less controversial, there are obvious questions regarding the neighborhood size and/or position: Would the results change if a larger/smaller diameter were used? What if the disk center was offset from the GA margin? While we leave the answers to these questions to future studies, we note here that, qualitatively, a 500 µm diameter neighborhood provides a reasonable balance between localization and averaging, and also represents a physiologically plausible region of influence.
Our GA growth model also has limitations. First, although our model has some basic biological plausibility, it surely represents an incomplete picture of the true spatiotemporal dynamics of GA growth. And, while the model generates local trajectories that are prima facia reasonable, the model, its parameters (i.e., α and β), and its resulting trajectories need to be calibrated and validated against ground truth measurements-which, considering the required frequency of imaging, and the eye-to-eye variability in GA growth patterns, are difficult to obtain. However, in general, we believe that our model captures overall trends in local GA growth rates. Moreover, it is reasonable to expect that our model most accurately generates GA growth trajectories for smaller GA growths-in cases with larger lesion growths, and particularly those involving complex merging of multiple foci, the array of plausible growth trajectories dramatically increases. Further development, calibration, and validation of our growth model promises to be an interesting direction for this and other applications. Another limitation of our GA growth modelling is that it does not consider GA growth rates in the context of their location and direction. Specifically, GA lesions exhibit, on average, different growth rates along different directions-for example, growing more slowly towards the fovea [51]. Therefore, comparing an X mm/year growth rate in the immediate vicinity of the fovea to an X mm/year growth rate far from the fovea may be misleading. Possible strategies to account for such variations are similar to those proposed to account for trends in CC impairment with distance from the fovea-namely normalizing growth rates based on a statistical growth atlas or grouping margin segments according to their position and/or growth direction.
In addition to the above limitations, there are more general limitations, shared, to varying degrees, by most other CC OCTA studies. Foremost, the observed correlations between local CC FD percentages and local GA growth rates are only biologically meaningful to the extent to which measured CC FD percentages and GA growth rates reflect the physiological CC flow impairments and GA growths rates. The adequacy of our measurements is complicated by physiological factors (e.g., vitreous opacities and ocular abberations), as well as by limitations and parameter choices in the imaging instrumentation, image acquisition, processing, and analysis. Exhaustive enumeration and description of these limitations and parameter choices is beyond the scope of this paper. However, to emphasize the complexity and potential for error, we provide an incomplete list of the steps involved in the analysis, along with their associated limitations and parameter choices: OCT imaging and parameters (e.g., incident power, A-scan rate, wavelength, optical spot size, vignetting, and motion artifacts); OCTA processing methods and parameters (e.g., normalization method, interscan time, number of repeated B-scans); OCTA analysis methods and parameters (e.g., segmentation artifacts, choice of projection range, OCTA binarization method and parameters); and GA tracing (e.g., slab used to trace GA margins, definition of GA). Variations in any of these processing steps, and/or poor choices of any of these parameters, can result in measurements that diverge from true biological values. To highlight these potential complications, we explore below two ways in which OCTA measurements of CC impairment can diverge from true physiological CC impairments.
One challenge in OCTA CC imaging is the coupling between the OCTA CC signal and the intensity of the backscattered light. For example, fewer backscattered photons will be collected from regions of the CC with overlying drusen than from regions of the CC without overlying drusen, which can cause differences in CC OCTA signals even if there are no differences in the CC blood flows. Typically, such effects are mitigated by normalizing the OCTA signal by the OCT signal using one of two approaches: normalized OCTA metrics or post-processing normalization. When a normalized OCTA metric is used, the outputted OCTA signal is, by construction, normalized by the OCT signal; this is the approach taken by the split-spectrum amplitude-decorrelation angiography (SSADA) algorithm [52], currently used by Optovue (Fremont, CA, USA), as well as the OCTA ratio analysis (OCTARA) algorithm [53], currently used by Topcon (Oakland, NJ, USA). In post-processing approaches, an unnormalized OCTA metric is used and OCTA images are normalized in a subsequent processing step, typically following en face projection. This approach, presented by Zhang et al. [28] and used in our present study, is applicable to OCTA data generated using the complex optical micro-angiography (OMAG c ) algorithm [6], which is currently used by Carl Zeiss Meditec (Dublin, CA, USA) instruments. In our view, both are reasonable approaches to decoupling the OCTA signal from the intensity of backscattered light. However, neither approach is perfect, and neither approach works well when the OCT signals are strongly attenuated, such as in the deeper choroid-in these regions there is simply insufficient backscattered light to reliably perform OCTA. A recent study by Ledesma-Gil et al. [54] suggests that the compensation of Zhang et al.-which has been also been used in other OCTA studies of the CC, such as those by Nassisi et al. [13] and Alagorie et al. [14]-can induce artifacts, for example by producing new, artifactual flow deficits, and that different CC slab positions change the appearance of these artifacts. While we cannot rule out that our findings were influenced by limitations in OCTA signal compensation, it is worth noting that our previous pilot study [11], which used an unnormalized and unthresholded analysis strategy [55], found no obvious, general correlations between local CC impairment and local GA growth rates. However, future studies using varied OCTA imaging and compensation techniques are needed.
Another example of a challenge arising in OCTA CC imaging is the dependency of OCTA signals in general [41,[56][57][58], and of OCTA CC FD deficits in particular [59], on the interscan time. In the present study, the 6.1 ms interscan time-typical for commercial instruments-likely does not capture more subtle CC blood flow impairments. However, as evidenced by the CC FD images of this study, the 6.1 ms interscan time does visualize many flow deficits around the margin. Given that these more severe flow deficits were not found to correlate well with local GA growth rates, it is not obvious why more subtle flow impairments would improve correlations. Moreover, our previous pilot study [11], which assessed both 1.5 ms and 3.0 ms interscan time OCTA, found no obvious, general correlations between local CC impairment and local GA growth rates. Nevertheless, future OCTA studies using different-and, in particular, shorter-interscan times are needed.

Conclusions
In this study we analyzed correlations between local CC FD percentages and local GA growth rates using SS-OCTA imaging, GA growth modeling, and spatial statistical methods. The null hypothesis (Pearson's ρ = 0) was rejected at a Benjamini-Hochberg false discovery rate of 0.2 in 15 of the 114 visit pairs, 11 of which exhibited positive correlations; even amongst these 11 visit pairs, correlations were modest (r in [0.30, 0.53]). Nevertheless, future studies using different instruments and quantification techniques are needed. The analysis framework presented in this study appears well suited to rigorously assessing correlations between local GA growth rates and other imaging biomarkers.