Classification of Visual Cortex Plasticity Phenotypes following Treatment for Amblyopia

Monocular deprivation (MD) during the critical period (CP) has enduring effects on visual acuity and the functioning of the visual cortex (V1). This experience-dependent plasticity has become a model for studying the mechanisms, especially glutamatergic and GABAergic receptors, that regulate amblyopia. Less is known, however, about treatment-induced changes to those receptors and if those changes differentiate treatments that support the recovery of acuity versus persistent acuity deficits. Here, we use an animal model to explore the effects of 3 visual treatments started during the CP (n = 24, 10 male and 14 female): binocular vision (BV) that promotes good acuity versus reverse occlusion (RO) and binocular deprivation (BD) that causes persistent acuity deficits. We measured the recovery of a collection of glutamatergic and GABAergic receptor subunits in the V1 and modeled recovery of kinetics for NMDAR and GABAAR. There was a complex pattern of protein changes that prompted us to develop an unbiased data-driven approach for these high-dimensional data analyses to identify plasticity features and construct plasticity phenotypes. Cluster analysis of the plasticity phenotypes suggests that BV supports adaptive plasticity while RO and BD promote a maladaptive pattern. The RO plasticity phenotype appeared more similar to adults with a high expression of GluA2, and the BD phenotypes were dominated by GABAAα1, highlighting that multiple plasticity phenotypes can underlie persistent poor acuity. After 2-4 days of BV, the plasticity phenotypes resembled normals, but only one feature, the GluN2A:GluA2 balance, returned to normal levels. Perhaps, balancing Hebbian (GluN2A) and homeostatic (GluA2) mechanisms is necessary for the recovery of vision.


Introduction
Since the earliest demonstrations that monocular deprivation (MD) during a critical period (CP) causes ocular dominance plasticity and acuity loss [1][2][3], this model has been used to deepen our understanding of the neural changes associated with amblyopia. There have been fewer studies, however, about cortical changes associated with the acuity deficits that often persist after treatment for amblyopia [4][5][6][7][8]. Here, we use an animal model to classify the expression patterns (phenotypes) of a collection of synaptic proteins that regulate experience-dependent plasticity and explored if treatments that promote good versus poor acuity reinstate CP-like plasticity phenotypes in the visual cortex (V1).
Many animal studies have highlighted the roles of glutamatergic and GABAergic mechanisms for regulating plasticity during the CP [9][10][11][12][13][14][15]. For example, the subunit composition of AMPA, NMDA, and GABA A receptors regulates the bidirectional nature of ocular dominance plasticity [16][17][18][19][20][21]. Some of the changes caused by MD include delaying the maturational shift to more GluN2A-containing NMDARs [22,23] and accelerating the expression of GABA A α1-containing GABA A Rs [20,23]. Together, those changes likely decrease signal efficacy and dysregulate the spike-timing-dependent plasticity that drives long-term depression (LTD) and weakens deprived-eye response [24]. Furthermore, silencing activity engages homeostatic mechanisms that scale the responsiveness of V1 neurons by inserting GluA2-containing AMPAR into the synapse [25]. Importantly, many of the receptor changes have been linked with specific acuity deficits [26,27] suggesting that visual outcomes may reflect changes to a collection of glutamatergic and GABAergic receptor subunits that together represent a plasticity phenotype for the V1.
Animal studies of amblyopia have also identified treatments that promote good versus poor recovery of acuity after MD. For example, reverse occlusion (RO) gives a competitive advantage to the deprived eye that promotes an ocular dominance shift, but the acuity recovered by the deprived eye is transient and can be lost within hours of introducing binocular vision [6][7][8]. Similarly, closing both eyes after MD to test a form of binocular deprivation therapy (BD) leads to poor acuity in both eyes that does not recover even after months of binocular vision [28]. In contrast, just opening the deprived eye to give binocular vision (BV) after MD appears to engage cooperative plasticity that promotes both physiological recovery [29] and long-lasting visual recovery in both eyes [27].
Here, we quantified the expression of glutamatergic and GABAergic receptor subunits in the V1 of animals reared with MD and then treated to promote either good visual recovery (BV) or persistent bilateral amblyopia (RO, BD). Next, we developed an unbiased high-dimensional analysis approach to identify plasticity features in the pattern of subunit expression and to construct plasticity phenotypes. Finally, we used cluster analysis to classify plasticity phenotypes associated with good versus poor acuity and analyzed those to determine which features suggest the recovery of adaptive versus maladaptive plasticity mechanisms.

Animals and Rearing
Conditions. All experimental procedures were approved by the McMaster University Animal Research Ethics Board. We quantified the expression of 7 glutamatergic and GABAergic synaptic proteins in the V1 of cats reared with MD from eye opening until 5 weeks of age and then given one of the 3 treatments: RO for 18 d, BD for 4 d, or BV for either short-term (ST-BV; 1 hr, 6 hrs) or long-term (LT-BV; 1 d, 2 d, or 4 d) (n = 7, 4 male and 3 female) ( Figure 1). The lengths of RO and BD were selected because they have well-documented and consistent visual changes that result in poor acuity in both eyes [7,8,30]. The BV periods were selected to match the lengths used previously to study rapid and dynamic changes caused by MD in both cat and mouse V1 [27,31,32]. Also, the shortand long-term BV groups were based on the data-driven analysis of protein expression described in detail below and that analysis placed the samples from ST-BV (1 hr or 6 hrs) versus LT-BV (1 d, 2 d, or 4 d) rearing conditions into separate clusters. The raw data collected previously [23] from animals reared with normal binocular vision until 2, 3, 4, 5, 6, 8, 12, 16, or 32 wks of age (n = 9 animals, 2 male and 7 female) or MD from eye opening (6-11 d) to 4, 5, 6, 9, or 32 wks (n = 8 animals, 4 male and 4 female) were used for comparison.
MD was started at the time of eye opening by suturing together the eyelid margins of one eye (5-0 Coated VICRYL Ethicon P-3) using surgical procedures described previously [8]. Sutures were inspected daily to ensure the eyelids remained closed. At 5 weeks of age, the period of MD was stopped and either BV was started by carefully parting the fused eyelid margins, RO was started by opening the closed eye and closing the open eye, or BD was started by closing the open eye. All of these surgical procedures were done using gaseous anesthesia (isoflurane, 1.5-5%, in oxygen) and aseptic surgical techniques.
At the end of the rearing condition, animals were euthanized using sodium pentobarbital injection (165 mg/kg, IV) and transcardially perfused with cold 0.1 M phosphatebuffered saline (PBS) (4°C; 80-100 ml/min) until the circulating fluid ran clear. The brain was removed from the skull and placed in cold PBS. A number of tissue samples (2 mm × 2 mm) were taken from the regions of the V1 representing the central (C), peripheral (P), and monocular (M) visual fields (Figure 1(c)). Each tissue sample was placed in a cold microcentrifuge tube, flash frozen on dry ice, and stored in a -80°C freezer.
Total protein concentrations were determined for each sample and a set of protein standards using the bicinchoninic acid (BCA) assay (Pierce, Rockford, IL, USA). A linear function was fit to the observed absorbance values of the protein standards relative to their expected protein concentrations. If the fit was less than R 2 = 0 99, the assay was redone. The slope and the offset of the linear function were used to determine the protein concentration of each sample, and then the samples were diluted to 1 μg/μl with sample (M260 Next Gel Sample loading buffer 4x, Amresco) and Laemmli buffer (Cayman Chemical). A control sample was made by combining a small amount from each sample to create an average sample that was run on every gel. Each sample was run twice in the experiment.
The appropriate secondary antibody conjugated to horseradish peroxidase (HRP) (1 : 2,000; Cedarlane Laboratories LTD, Hornby, ON) was applied to membranes for 1 hour at room temperature, then blots were washed in PBS (3 × 10 min). Bands were visualized using enhanced chemiluminescence (Amersham, Pharmacia Biotech, Piscataway, NJ) and exposed to autoradiographic film (X-Omat, Kodak, Rochester, NY). After each exposure, blots were stripped (Blot Restore Membrane Rejuvenation kit (Chemicon International, Temecula, CA, USA)) and probed with the next antibody so each blot was probed for all 7 antibodies (Figure 1(d)).

Analysis of Protein
Expression. The autoradiographic film and an optical density wedge (Oriel Corporation, Baltimore, MD) were scanned (16 bit, AGFA Arcus II, Agfa, Germany), and the bands were identified based on molecular weight. The bands were quantified using densitometry, and the integrated grey level of the band was converted into optical density units (OD) using custom software (MATLAB, The MathWorks Inc., Natick, Massachusetts). The background density between the lanes was subtracted from each band, and the density of each sample was normalized relative to the control sample run on each gel (sample band density/control band density).  Figure 1: Study design diagram. Timelines for the rearing conditions used in this study. (a) Normal visual experience and monocular deprivation (MD), (b) treatment conditions (RO, BD, and BV) after MD to 5 wks. Filled bars indicate that an eye was closed. Black arrows indicate the age of animals used in the study. A timeline for the critical period (CP) in cat visual cortex [34] highlights the peak of the CP between 4 and 6 weeks of age. (c) Map of the sampling regions in the V1 representing the central (C, n = 2), peripheral (P, n = 8), and monocular (M, n = 2) visual fields. (d) Representative bands from the Western blots for the 7 proteins quantified in the study and the molecular weights (kDa).
The data were normalized relative to the average expression of the 5 wk normal cases. Table 2 summarizes the number of tissue samples and replication of runs for the 5 wk normal, 5 wk MD, and recovery conditions across the 3 regions of the V1 and 7 proteins that were studied. Descriptions of the expression for the individual proteins in each of the conditions can be found in [35]. Those univariate comparisons confirmed the complex nature of these data and led us to develop and implement the data analysis workflow that is summarized in Figure 2. 2.5. Protein Network Analysis. A network analysis of protein expression was done for each rearing condition by calculating the pairwise Pearson's R correlations among the 7 proteins using the rcorr function in the Hmisc package in R [36]. The networks were visualized as correlation matrices (heatmap2 function in gplots [37]), and the proteins were ordered using the dendextend [38] and seriation [39] packages to place proteins with similar patterns of correlations nearby in the dendrogram. Significant correlations were identified using the Bonferroni-corrected p values and indicated by asterisks on the cell in the correlation matrix.
2.6. Principal Component Analysis. We used principal component analysis (PCA) to reduce the dimensionality of the data, identify potential biological features, and create plasticity phenotypes. We applied the PCA following the procedures we used previously [23,40,41] and included data from all of the normal animals and MDs as well as the 3 recovery conditions. We assembled the protein expression for GluA2, GluN1, GluN2A, GluN2B, GABA A α1, GABA A α3, and synapsin into an mxn matrix. The m columns  represented the 7 proteins, and the n rows were the average protein expression for each of the 12-14 samples from an animal. For a few of the rows, data was missing from a single cell, and so those samples were omitted for a total of n = 279 rows in the matrix and 1,953 observations. The data were centered by subtracting the mean column vector and applying singular value decomposition (SVD) to calculate the principal components (RStudio). SVD represents the expression of all 7 proteins within a single tissue sample as a vector in a high dimensional space, and the PCA identifies variance captured by each dimension in that "protein expression space." The first 3 dimensions accounted for 82% of the total variance and were used for the next analyses.
We plotted the basis vectors for the first 3 dimensions (Dim) and used the weight, quality (cos 2 ), and directionality of each protein, as well as known protein interactions, to help identify potential biological features accounting for the variance. We identified 9 potential features, calculated those features for each sample, and correlated each feature with Dim1, Dim2, and Dim3 to create a correlation matrix (see results). The p values for the correlations were Bonferroni corrected, and significant correlations were used to identify the features that would be part of the plasticity phenotype.
Eight of the features were significantly correlated with at least one of the first 3 dimensions. A measure associated with the E:I balance was not significantly correlated with the dimensions, and so it was not included in the tSNE or cluster analysis. The E:I measure, however, was used for analyzing the composition of the clusters and as a component of the plasticity phenotype because of the importance of the E:I balance for experience-dependent plasticity.

tSNE Dimension
Reduction and Cluster Analysis. The average expression for the 8 features (Table 3) was compiled into an mxn matrix, with m columns (m = 8) representing the significant features and n rows representing each sample from the 3 V1 regions (central, peripheral, and monocular) for 5 wk normal, 5 wk MD, RO, BD, and BV animals (n = 109). t-distributed stochastic neighbor embedding (tSNE) was used to reduce this matrix to 2dimensions (2D). tSNE was implemented in R [42], and the tSNE output was sorted using k-means to assign each sample to a cluster. To determine the optimal number of clusters (k), we calculated the within-groups sum of squares for increasing values of k, fit a single-exponential tau decay function to those data, found the "elbow point" at 4τ which was 6, and used that as the optimal number of clusters. The clusters were visualized by color-coding the dots in the tSNE plot, and the composition of the clusters was analyzed.
To facilitate analysis of the tSNE clusters, we grouped the BV cases into short-term BV (1 hr and 6 hr) (ST-BV) or long-term BV (1 d, 2 d, and 4 d) (LT-BV), color-coded the samples by rearing condition, and used different symbols to indicate the V1 region. For each cluster, we annotated the composition based on the rearing condition of the samples to create "subclusters" (e.g., LT-BV 1) that were used for the next analyses.
We evaluated the similarity/dissimilarity among the subclusters by calculating the pairwise correlations (Pearson's R)  between subclusters using the features identified by the PCA as input to the R package rcorr. The correlations were visualized in a matrix with the cells color-coded to indicate the strength of the correlation [37]. The order of the subcluster in the matrix was optimized using hierarchical clustering, and a dendrogram was created based on the pattern of correlations (using dendextend and seriation packages in R) so that subclusters with strong correlations were nearby in the dendrogram.

Visualization and Comparison of Plasticity Phenotype.
The features identified in the PCA were used to indicate the plasticity phenotype for each of the subclusters. In addition to the 8 significant features, the E:I measure was included in the visualization of the plasticity phenotype. The features were color-coded using grey scale for the 3 protein sum features and a color gradient (red = -1, yellow = 0, and green = +1) for the 6 protein indices. The plasticity phenotypes were displayed as a stack of color-coded bars with one bar for each feature. For the subclusters, the plasticity phenotypes were ordered by the dendrogram to facilitate comparison among subclusters that were similar versus dissimilar. We also calculated the plasticity phenotypes for the full complement of normally reared and MD animals and displayed those in a developmental sequence to facilitate age-related comparisons with the recovery subclusters. Finally, we did a bootstrap analysis to determine which features of the plasticity phenotypes were different from 5 wk normals and used Bonferroni correction to adjust the significance for the multiple comparisons. This analysis was displayed in 2 ways: first, each of the 9 feature bands for the dendrogram-ordered subclusters was color-coded with white if it was not different, red if it was greater, and blue if it was less than 5 wk normals; second, boxplots were made to show the value for each of the 9 features and to identify the subclusters that were different from 5 wk normals.
A detailed description of the network analysis, PCA, tSNE, clustering, and phenotype construction, along with the example code for each of these steps, can be found in [43].

Modeling Population Receptor Decay Kinetics for NMDARs
and GABA A Rs. The subunit composition of NMDARs and GABA A Rs determines the decay kinetics of the receptor [44,45], and so we used that information to build a model for the decay kinetics of a population of receptors for each of the rearing conditions. The decay kinetics of the most common NMDAR composition, triheteromeric receptors containing GluN2A and 2B, is 50 ms ± 3 ms, while diheteromers NMDARs containing only GluN2B are slower (2B = 333 ms ± 17 ms) and those containing only GluN2A are faster (2A = 36 ms ± 1 ms) [44]. The decay kinetics of GABA A Rs with both α1 and α3 subunits is 49 ms ± 23 ms while receptors with only the α3 subunit are slower (129 0 ms ± 54 0 ms) and only α1 are faster (42 2 ms ± 20 5 ms) [45].
We used the relative amounts of GluN2A and 2B, or GABA A α1 and α3, as inputs to the model. Receptors containing GluN2A and 2B or GABA A α1 and α3 are the most common in the cortex, so the model maximized the number of these pairs which was limited by the subunit with less expression. The remaining proportion of the highly expressed subunit was divided by 2 and used to model the number of pairs for those receptors (2A:2A or 2B:2B; α1:α1 or α3:α3) in the population. The population decay kinetics were then modeled by inserting the relative amounts of the subunits into these formulas: Table 3: Formulas and Pearson's R correlation between the features and principal components. The formulas for PCA-identified features, including protein sums ( Figure 5) and receptor indices ( Figure 6), along with corresponding correlation (R 2 ) values for each of the first 3 principal components. The GluN1:GluA2 and GABA A R sum:GlutR sum were not significantly correlated with any of these 3 components.
First, we plotted scattergrams of the average NMDAR and GABA A R decay kinetics for normal animals and each treatment condition. The development of decay kinetics for normal animals was described using an exponential decay function, while changes in kinetics with increasing lengths of BV were fit by either exponential decay or sigmoidal curves. Then, we compared the relationship between NMDAR and GABA A R kinetics by plotting both on one graph.

Statistical Analyses.
We used the bootstrap resampling method to compare the features because it is a conservative approach to analyzing small sample sizes when standard parametric or nonparametric statistical tests are not appropriate. Here, bootstrapping was used to estimate the confidence intervals (CI) for each feature in the subcluster, and a Monte Carlo simulation was run to determine if the 5 wk normal subcluster fell outside those CIs. The statistical software package R was used to simulate normal distributions with 1,000,000 points using the mean and standard deviation from the subcluster. Next, a Monte Carlo simulation was randomly sampled with replacement from the simulated distribution n times, where n was the number of observations made from the normal subcluster. The resampling procedure was repeated 100,000 times to determine the 95%, 99%, and 99.9% CIs. The subcluster feature was considered significantly different from normal (e.g., p < 0 05, p < 0 01, or p < 0 001) if the feature mean fell outside these CIs. When a subcluster was significantly greater than normal (p < 0 05), the boxplot was colored red; when it was less than normal (p < 0 05), the boxplot was colored blue; and if it was not different from normal (p > 0 05), the boxplot was colored grey.
All of the bootstrap statistical comparisons for the plasticity features (Table 5-1 and 6-1) are presented in the Supplemental material.
The p values for Pearson's correlations were calculated using the rcorr package [36], and the significance levels were adjusted using the Bonferroni correction for multiple comparisons. Pearson's Rs and p values for the protein networks (Table 3-1), plasticity features with PCA dimensions (Table 4-1), and association between clusters (Tables 8-1 We tested if recovery during BV followed either an exponential decay or a sigmoidal pattern by fitting curves to the data using Kaleidagraph (Synergy Software, Reading, PA). Significant curve fits were plotted on the graphs to describe the trajectory of recovery.

Analyzing the Pairwise Similarity between Protein
Expression Profiles. First, we wanted to identify pairs of proteins with similar or opposing expression profiles and compare them among the rearing conditions. For each condition, we collapsed the data from the 3 regions of the V1, calculated the matrix of pairwise correlations between the 7 proteins, ordered the protein correlations using a hierarchical dendrogram, and used 2D heatmaps to visualize the correlations ( Figure 3). The order of proteins in the dendrogram indicated which ones had similar (e.g., on the same branch of the dendrogram) or different patterns of expression, and the color of the cell illustrated the strength of the correlation. For 5 wk normal animals (Figure 3(a)), there were strong positive correlations (red cells) among all proteins except GluN2A, which was weakly correlated and not clustered with the other proteins. A different pattern of correlations was found after MD (Figure 3(b)); here, glutamatergic proteins were weakly or even negatively correlated (blue cells) with GABA A α1, GABA A α3, and synapsin. These results suggest that MD drives a decoupling of these excitatory and inhibitory mechanisms. RO also separated glutamatergic and GABAergic proteins into different clusters at the first branch (Figure 3(c)); however, the correlations were weaker, suggesting that RO reduced the MD-driven decoupling of these mechanisms. After BD, the correlation matrix had mostly positive correlations (Figure 3(d)) except for synapsin which was negatively correlated and not clustered with the other proteins. BV treatment highlighted the dynamic nature of this recovery (Figures 3(e)-3(i)). Just 1 hr of BV was enough to change the correlation matrix from the MD pattern, but even after 4 d of BV, the correlation matrix still appeared different from the normal 5 wk pattern of correlations.
These matrices suggest different patterns of correlations depending on the condition, but this analysis treats each comparison with the same weighting and it is likely that some proteins contribute more than others to the variance in the data. To assess this, we used the PCA to identify individual proteins and combinations of proteins that capture the variance in the data and may represent plasticity features reflecting differences among the treatment conditions.

Using Principal Component Analysis to Reduce
Dimensionality and Identify Plasticity Features. We used the PCA to reduce the dimensionality, transform the data, and find features that define the covariance among the proteins. An mxn matrix was made using protein expression, where the m columns were the 7 proteins and the n rows (109) were the tissue samples from all the animals and regions of the V1 used in this study. This matrix was analyzed using singular value decomposition (SVD), and the first 3 dimensions explained most of the variance (82%) in the data (Dim1: 54%, Dim2: 18%, and Dim3: 10%) (Figure 4(a)).
To understand which proteins contributed to each dimension, we addressed the quality of the representation for each protein using the cos 2 metric and found that the glutamatergic proteins were well represented by Dim1, GABA A α1 by Dim2, and GluA2 and GluN2B by Dim3, but synapsin and GABA A α3 were weakly represented in the first 3 dimensions (Figures 4(b) and 4(c)). Next, we compared the vectors for each protein (Figures 4(d) and 4(f)) and the PCA space occupied by the rearing conditions (Figures 4(e) and 4(g)). The protein vectors show that GluN1, GluN2A, GluN2B, and GluA2 extended along Dim1, GABA A α1 along Dim2, and GluA2 and GluN2B were in different directions along Dim3. The PCA space occupied by the conditions suggest some differences: BD was separated on Dim2 in the same direction as GABA A α1, but the center of gravity for the other conditions overlapped the space occupied by normal samples.
The overlap among conditions raised the possibility that higher dimensions may separate the conditions. To begin to assess higher dimensional contributions, we examined the basis vectors (Figure 4(h)) and the correlations between individual proteins and PCA dimensions (Figure 4(i)) to identify combinations of proteins that might reflect higher dimension features. For example, all proteins had positive amplitudes for the Dim1 basis vector (Figure 4(h)), and positive correlations with Dim1 (Figure 4(i)) suggested that protein sums may be higher dimensional features. In addition, weights for GluN2A and GABA A α1 on Dim2 were opposite, suggesting that when one protein increased the other decreased, and this could be a novel feature of these data. Continuing with this approach, we identified 9 putative plasticity features: protein sums (all protein sum, GlutR sum, and GABA A R sum) or indices (GlutR:GABA A R, Glu-N2A:GluN2B, GABA A α1:GABA A α3, GluN2A:GABA A α1, GluA2:GluN2B, and GluN2A:GluA2). All of the protein sums and 4 of the indices were features not analyzed with the univariate statistics; however, each had a strong biological basis in previous research. For example, the new indices paired the mature GluN2A with the mature GABA A α1 or GluA2 subunit and GluN2B with GluA2 which is known to regulate the development of AMPARs [46]. Finally, we calculated the 9 features and determined if at least one of the first 3 dimensions was correlated with the features (Figure 4(j)). Only the GlutR:GABA A R balance was not correlated with any of the first 3 dimensions, but because those mechanisms are related to the E:I balance [47], we included that measure in the next analysis.

Comparing Plasticity Features.
We plotted the plasticity features and saw that the GlutR and GABA A R sums and indices identified various differences among the treatment conditions ( Figures 5 and 6). There were, however, consistent changes after BV in the binocular regions with a loss of the total amount of GABA A R expression (44% ± 12) and a shift of the GlutR:GABA A R balance to favor GlutR ( Figure 5(d)). The remaining indices in the feature list also identified differences ( Figure 6) including the GABA A α1:GluN2A balance shifting to more GluN2A after BV (in binocular regions) but more GABA A α1 after BD. RO flipped the 2A:2B balance GluA2  Six clusters were visualized with tSNE (Figure 7), and the composition of the clusters was analyzed to determine the V1 regions and rearing conditions in each cluster. Cluster 1 was the largest with 39 samples (C = 26%; P = 54%; and M = 21%) and had the greatest number of samples from the central region (Figures 7(b) and 7(d)). Cluster 3 also had samples from the central, peripheral, and monocular regions while clusters 4, 5, and 6 were dominated by peripheral samples with few or no central region samples. Thus, there was some clustering by the V1 region, but more apparent clustering emerged when the samples were color-coded by rearing condition (Figures 7(c) and 7(d)). All but one of the normal samples were in cluster 1, all of the RO samples were in cluster 2, most of the BD samples were in cluster 3 with a few in cluster 6, and most of the MD samples were in clusters 1 or 3. The BV samples, however, were found in 5 of the clusters with the greatest number of BV central samples (83%) grouped with normals in cluster 1. GluN2B GluN1 GluA2   Further analysis of cluster 1 showed that the majority of LT-BV and ST-BV samples from the central region clustered with the normals (Figure 7(d)). Interestingly, some of the MD samples were also in cluster 1; however, those samples were from the peripheral and monocular regions which are known to be less affected by MD than the central region [48]. Together, these results show that the data are clustered and that the clustering was driven by both the rearing condition and the region of the V1.

Correlating Plasticity Features among Subclusters.
We annotated the samples in each cluster using the rearing condition and V1 region and used that information to identify 13 subclusters where at least one region per condition had n ≥ 2 and >20% of the samples in that cluster (Figure 7(d), black font). A correlation matrix was calculated (Figure 8) to assess the similarity between subclusters (see Supplemental Table 8

Constructing Plasticity Phenotypes and Comparing among the Subclusters.
To compare the pattern of plasticity features among the subclusters, we visualized the average for each feature as a color-coded horizontal band, stacked the bands to illustrate the pattern that we called the plasticity phenotype (Figure 9(a)), and ordered the phenotypes using the same dendrogram as the correlation matrix (Figure 9(b)). In addition, we visualized the plasticity phenotypes for normal development and MD (using the data from [23]) to compare the treatment subclusters with a broad range of ages that had developed with either normal or abnormal visual experience (Figures 9(c) and 9(d)).
Inspection of the plasticity phenotypes identified some obvious and other subtler differences among the subclusters (Figure 9(b)). Indeed, the pattern of red and green bands in the BD phenotypes was different from 5 wk nor-mals ( Figure 9) and showed the shift to more GABA A α1 and less GluN2A. For the RO subcluster, the light grey bands and number of green bands identified loss of protein expression and a shift to more GluN2A than 2B and more GluA2 than 5 wk normals. The RO pattern, however, appeared similar to an older (e.g., 12 wk) normally reared animal suggesting that RO may accelerate maturation of these proteins. Thus, these BD and RO treatments led to distinct plasticity phenotypes.
The pattern of red and green bands in the plasticity phenotype for LT-BV and some of the ST-BV subclusters (ST-BV1, ST-BV5) appeared similar to the 5 wk normals (Figure 9(b)), but many of the features were still significantly different from the age-matched normals (Figure 10(a), Supplemental Table 10-1). Nonetheless, these subclusters had some consistent differences with less GABA A Rs and more GluN2B than 5 wk normals. Interestingly, one of the novel features found by the PCA, the GluN2A:GluA2 balance, was the only measure where all of the LT-BV subclusters were not different from 5 wk normals, but both   (Figure 4(j)) and plotted for each region of the V1 (columns). The conventions are the same as in Figure 5. For exact p values, Pearson's R, and equations for the curve-fits, see Supplemental Table 6-1. RO and BD were different. Thus, this visualization of the plasticity phenotypes illustrated that the pattern promoted by BV, and LT-BV in particular, was most similar to the normal CP phenotype.

Modeling NMDAR and GABA A R Population Kinetics.
The subunit composition of NMDARs and GABA A Rs helps to regulate the threshold for experience-dependent plasticity, in part by controlling receptor kinetics [44,45]. We used the  information about receptor kinetics with different subunit compositions to make a model that predicts the average population kinetics and applied it to normal development and the rearing conditions studied here. First, we transformed the 2A:2B and α1:α3 balances into predicted population kinetics (see Methods) and plotted the normal postnatal development (Figures 11(a) and 11(b)). Both the NMDA and GABA A kinetics rapidly speed up between 2 and 6 weeks of age. Next, we compared the predicted kinetics among the rearing conditions (Figures 11(c) and 11(d)). The pattern of results is necessarily similar to the balances presented for the indices ( Figure 6); however, the predicted kinetics suggests a compression of differences between conditions when the balances favor the mature subunits (2A or α1) versus an accentuation of differences with much slower kinetics when the immature subunits (2B or α3) dominated.
To address how treatment-induced changes to NMDAR and GABA A R composition might affect the relationship between glutamatergic and GABAergic transmission timing, we made XY scatterplots using the predicted kinetics ( Figure 11(e)). During normal development (black line), both balances progressed from slow kinetics at 2 wks to faster kinetics through the peak of the CP (Figure 11(e), yellow zone; 4-6 wks) to reach adult levels. The NMDAR:GABA A R kinetics for MD, RO, and BD fell outside the window predicted for the normal CP but in different directions. MD had slower NMDAR (C: 135 ms ± 16 ms; P: 121 ms ± 12 ms; and M: 146 ms ± 27 ms) but faster GABA A R kinetics (C: 47 ms ± 0 3 ms; P: 48 ms ± 1 ms; and M: 51 ms ± 4 ms), RO had faster NMDAR (C: 46 ms ± 0 8 ms; P: 46 ms ± 0 4 ms; and M: 46 ms ± 0 2 ms) but normal CP range for GABA A R (C:54 ms ± 6 ms; P:51 ms ± 2 ms; and M: 48 ms ± 0 2 ms), and BD had faster GABA A R (C: 46 ms ± 0 9 ms; P: 44 ms ± 0 2 ms; and M 45 ms ± 0 2 ms) but normal CP range NMDAR kinetics in the central region only (C: 61 ms ± 12 ms; P: 130 ms ± 12 ms; and M: 155 ms ± 27 ms).
The introduction of BV caused a progressive change in the predicted NMDAR:GABA A R kinetics suggesting an initial speeding up of the NMDAR kinetics over the first 1 d to 2 d followed by a slowing of the GABA A R kinetics, especially in the central region. Taken together, the predicted NMDAR:GABA A R kinetics provided additional evidence  Figure 7(d) and annotated here using the rearing condition, cluster (1-6), and region of the V1. Hierarchical clustering was used to order the data so that subclusters with strong correlations were nearby in the matrix. The subclusters formed 5 groups using the height of the dendrogram that is denoted by a change in the color of the dendrogram. The dotted black line in the dendrogram highlights the path to the normal subcluster. The black lines in the matrix identify the 5 groupings of the subclusters. For exact Bonferroni-corrected p values, see Supplemental that BV shifts protein expression towards a normal CP balance, but none of the treatments reinstated normal kinetics.

Discussion
Here, we studied a set of glutamatergic and GABAergic receptor subunits in the V1 that regulate plasticity and explored classifying treatments that cause either persistent bilateral amblyopia (RO or BD) or good acuity in both eyes (BV). Not surprisingly, there was a complex pattern of changes that varied by treatment and region within the V1. Applying a new analysis approach, however, using the PCA and cluster analysis, identified higher dimensional features and subclusters with different plasticity phenotypes for GluN2B -GluN2A GluN2B -GluA2 Adult 16wk 12wk 8wk 6wk 5wk 4wk 3wk 2wk

Protein sum
GlutR sum GABA A R sum GABA A R sum -GlutR sum GABA A 훼1 -GluN2A GluN2B -GluN2A GluN2B -GluA2 GluN2B -GluA2 Monocular deprivation 4wk 5wk 6wk 9wk Adult Figure 9: Visualizing the plasticity features and phenotypes for each subclusters. (a) We visualized the plasticity features as a stack of colorcoded horizontal bars that together comprise the plasticity phenotype. The 3 grey scale bars represent the protein sums, and the 6 red-green color-coded bars represent the protein indices identified by the PCA. (b) The plasticity phenotypes were calculated for each subcluster and ordered using the same dendrogram as described in Figure 8. (c) For comparison, the plasticity phenotypes were calculated using previously published data [23] for normal development (2 -32 wks) (d) and animals MDed from eye open until either 4, 5, 6, 9, or 32 wks.
GluN2B -GluN2A GluN2B -GluA2 GluN2A -GluA2 Cortical region treatments that promote good versus poor recovery of acuity. The LT-BV plasticity phenotypes were closest to the normal CP pattern while the RO phenotype appeared more similar to an older pattern dominated by GluA2. In contrast, the BD phenotypes were dominated by GABA A α1 making it distinct from RO and illustrating that multiple plasticity phenotypes can underlie persistent bilateral amblyopia. The PCA identified an understudied feature, the balance between mature glutamate receptor subunits (GluN2A:GluA2 balance), as a marker that might differentiate treatments supporting good acuity (BV), from those that lead to persistent bilateral amblyopia (RO, BD). Finally, modeling kinetics for NMDAR and GABA A R provided additional evidence that BV can return CP-like balances, especially in the central region of the V1.

Study Limitations and Design.
The exploratory nature of the design used here was limited because the small number of animals used leaves unanswered how much variation there is in response to the treatments. The visual manipulations (MD, RO, BD, and BV), however, are known to cause consistent changes in visual perception [7,8,[49][50][51][52], physiology [7,29,31,53], and molecular mechanisms [23,27,[54][55][56][57][58][59] that have been reliably measured in a number of laboratories using the cat to study visual system plasticity. Thus, these treatment-induced changes provide an understanding about the pattern of recovery that will be useful for formulating new hypotheses about the links between these proteins and persistent amblyopia. The study design had some strengths including that (i) the animal model has excellent spatial vision, with a central visual field, so we could compare changes in the regions of the V1 that represent different parts of the visual field [27], (ii) the treatments were initiated and completed within the CP [34], (iii) there is detailed information about the recovery of physiology for RO and BV [7,29,32,53] and acuity for all 3 treatments [7,8,27,29,30], (iv) both RO and BD cause persistent bilateral amblyopia [8,30], and (v) the treatments engage different forms of experience-dependent plasticity (RO: competitive; BD: cooperative with degraded visual input; and BV: cooperative with normal visual input).
We observed that only one feature (GluN2A:GluA2 balance) returned to normal after LT-BV treatment raising the hypothesis that it is necessary for good recovery. We were not able to test that question because the molecular tools are not available for manipulating proteins in the cat cortex so it will be important to replicate that finding in the mouse and then test the question by directly manipulating those proteins. In addition, a large number of other treatments have been tested to improve recovery after MD, including a brief period of darkrearing [30,60], fluoxetine administration [61], perceptual learning [27,62], or targeting specific molecular mechanisms (e.g., perineuronal nets [63]). Undoubtedly, the timing, length, and type of treatment influence recovery, but the conditions used here were necessarily limited because of the labor-intensive nature of this study. Notwithstanding these limitations, the plasticity phenotypes identified RO and BD as different from each other and from normals, but the LT-BV subclusters were remarkably similar to the 5 wk normal pattern.
Finally, the design took advantage of the reliability and multiplexing capabilities of Western blotting to obtain a large dataset of plasticity proteins from multiple V1 regions and rearing conditions. Western blotting, however, does not provide information about the cell types, layers, cortical columns, or subcellular localization of these proteins that would reveal which circuits are involved in recovery or persistent amblyopia. Even without that information, the application of high dimensional analyses led to the characterization of features and treatment-based clusters with unique plasticity phenotypes. The phenotyping approach developed here is scalable for studying more proteins or genes, cortical areas, and treatment conditions. Taken together, we think that  this approach can be used in other animal models where molecular tools can be combined with visual testing to identify the features and phenotypes necessary for optimal visual recovery.

BV Promoted Recovery of CP-Like Plasticity Phenotype and
Identified GluA2:GluN2A as a Balance That Differentiated BV Treatment. We explored BV treatment because it pro-motes long-lasting recovery of good acuity in both eyes [27], and those findings are similar to promising results of binocular therapies for treating amblyopia in children [64]. Furthermore, there is good physiological recovery with BV [29,32]. Thus, it was not surprising to find that LT-BV subclusters had the strongest correlations with normals or that those subclusters had CP-like phenotypes. However, only one of the features, the GluA2:GluN2A balance, returned     to normal levels. Those findings suggest that it may not be necessary to recapitulate every detail of the normal phenotype to support good visual recovery and that the GluA2:-GluN2A balance may be a characteristic feature for tracking functional recovery. Although that balance is not commonly quantified, both proteins are critical components of mechanisms regulating experience-dependent plasticity, and that balance might signify the adaptive engagement of multiple plasticity mechanisms. For example, the delayed increase in visual responses during ocular dominance plasticity is part of a homeostatic plasticity mechanism regulated by trafficking GluA2-containing AMPARs into the synapse [65,66]. Meanwhile, the initiation of ocular dominance plasticity requires GluN2A expression [22], and when GluN2A is deleted or reduced, MD does not depress deprived eye responses but instead causes enhancement of activity driven by the open eye [21]. Our finding that LT-BV returned a CP-like GluA2:GluN2A balance suggests that BV may prime GluN2A-dependent Hebbian plasticity to consolidate deprived-eye connections while GluA2-dependent homeostatic plasticity enhances deprived-eye responsiveness without triggering runaway excitation [67][68][69][70][71]. Thus, the GluA2:GluN2A balance may reflect the idea that during BV treatment the nondeprived eye acts as a teacher guiding both cooperative and competitive plasticity mechanisms [29].

RO versus BD Plasticity Phenotypes.
Because RO and BD treatments cause persistent bilateral amblyopia [7,8,30], we expected these conditions to have abnormal phenotypes. We were surprised, however, to find very different phenotypes for these conditions, showing that more than one plasticity phenotype can account for persistent acuity deficits. RO samples were in a single cluster dominated by an overabundance of GluA2 and more GluN2A than 2B. Together, those changes made the RO phenotype appear more similar to an adult than the CP pattern. The increase in GluA2 was in sharp contrast to the loss after BV treatment and suggests that RO may scale up AMPAR-dependent homeostatic mechanisms to drive recovery [25] without engaging NMDAR-dependent mechanisms to consolidate those changes [72]. Since AMPAR-mediated homeostasis promotes rapid but transient gains in responsiveness [25,65,[73][74][75][76], this might explain the labile acuity recovered with RO [7,8]. Interestingly, the overrepresentation of GluA2 promoted by RO implicates the dense expression of GluA2-containing synapses at feedback connections onto parvalbumin-positive (PV+) neurons [77]. The feedforward connections onto PV+ neurons may also be involved in RO circuit abnormalities because the labile acuity and early shift to GluN2A after RO are similar to changes found in MeCP2 KOs where an abnormally early shift to GluN2A at synapses onto PV+ neurons that halts acuity development [78,79]. Taken together, these findings provide preliminary evidence that RO may leave behind feedforward (GluN2A subunits) and feedback abnormalities (GluA2) in PV+ neuron circuits in the V1.
Although various models of neural plasticity predict that decreasing firing rates will enhance plasticity, that idea has not translated to using BD treatment to improve recovery from MD [30]. BD for weeks or months during the CP has a range of effects on the V1 including enhancing the appearance of cytochrome oxidase blobs [80], weakening stimulus-evoked activity of PV+ neurons [81], and delaying the developmental increase in the GAD65 expression [82]. Here, we found that a few days of BD treatment caused an abnormal increase in the expression of GABA A α1 throughout the V1 and a shift to more GluN2A in the central region. GABA A α1 receptors are found on pyramidal cell bodies where PV+ neurons synapse, and they serve as regulators of ocular dominance plasticity [20] and the window for coincident spike-time-dependent plasticity [24]. A recent study has shown that the loss of PV+ activity caused by BD depends on GABA A α1 mechanisms and that blocking this subunit increases BD-evoked activity allowing for LTP of PV+ neurons [83]. Our observation of increased GABA A α1 expression suggests that BD treatment may further reduce visually evoked activity in the V1 that is compounded by the shift to more GluN2A reducing the availability of the NMDAdependent mechanism needed to consolidate visual recovery.

Modeling Recovery of NMDAR and GABA A R Kinetics.
Our modeling of population kinetics suggests that different physiological changes accompany the 3 treatments. During normal development, the increases in NMDAR and GABA A R kinetics progress in concert. Physiological studies [84] and our modeling show that this fine balance is decoupled by MD because the delayed shift to GluN2A has slower NMDAR kinetics, but the premature increase of GABA A α1 has faster GABA A R kinetics. Neither RO nor BD treatment corrected that decoupling and the modeling suggests that those treatments accelerate the shift to faster adult-like kinetics for NMDARs after RO or GABA A Rs after BD. Modeling the kinetics for BV treatment identified 2 phases of recovery especially in the binocular regions of the V1. First, between 0 and 2 days of BV, there was a rapid increase in the predicted NMDAR kinetics that was similar to changes that occur between 2 and 4 weeks of age in normal cats. Second, between 2 and 4 days of BV, there was a slowing of the predicted GABA A R kinetics and that was opposite to the normal developmental increase in kinetics. These sequential phases of BV treatment do not recapitulate normal development. These results raise the question of whether the BV-driven increase in NMDAR kinetics needs to reach a certain level before triggering the slowing of GABA A R kinetics to rebalance these mechanisms. This modeling, however, was based on population data about the expression of the receptor subunits and cannot be extrapolated to individual receptors. Nonetheless, the rapid changes with BV treatment suggest that some aspects of normal development may be missed, and it will be important to determine what those are.

How Might These Plasticity Phenotypes be Used for Developing and Testing Treatments for Persistent Amblyopia?
The distinct plasticity phenotypes classified for RO and BD treatments provide preliminary evidence that multiple neural changes can account for persistent amblyopia and highlight the need to know which mechanisms to target when trying to engage neuroplasticity mechanisms to improve acuity. Whether the treatment should focus on AMPARs, NMDARs, GABA A Rs, or some combination of those receptors will depend on the underlying plasticity phenotype. Insights into those questions can be addressed in animal models using modern molecular tools and vision testing, but translating those findings into treatments for humans will depend on noninvasive ways to determine an individual's plasticity phenotype. For example, magnetic resonance spectroscopy has been used to measure changes in glutamate or GABA concentrations in human V1 after different types of visual experience (e.g., MD [85]), and receptor expression can be quantified by radioligands labeled for SPECT and PET [86]. New molecular imaging techniques hold the promise of even greater detail with the ability to measure the concentration of receptor subunits [87][88][89]. That information may be comparable to protein analysis in animal models and suitable for constructing plasticity phenotypes for human V1 to facilitate the translation of new treatments. Furthermore, behavioral paradigms linked with specific plasticity mechanisms (e.g., stimulus-selective response plasticity [90]) may further aid in characterizing human plasticity phenotypes. Thus, selecting a treatment to prevent or correct persistent amblyopia may benefit from in vivo steps to classify an individual's plasticity phenotype.

Conclusions
This exploration of glutamatergic and GABAergic receptor subunit changes in the V1 after treatment that promotes either good (BV) or poor (RO, BD) recovery of vision provides a better understanding of the complexity of this problem. Of the treatments studied here, only BV provided evidence for recovery of a CP-like plasticity phenotype in the V1. However, only one feature, the GluA2:GluN2A balance, returned to normal levels after BV, and that balance is noteworthy because the proteins are regulators of homeostatic and Hebbian plasticity, respectively. The modeling of NMDAR and GABA A R kinetics suggests two stages for BV recovery: a rapid increase in NMDAR kinetics, followed by slowing of the predicted GABA A R kinetics which together move that balance into the CP range. We identified features of the plasticity phenotypes that may guide future studies on persistent amblyopia to look for high levels of GluA2 and GluN2A following RO and high levels of GABA A α1 after BD treatment. Finally, the plasticity phenotyping is a good approach for uncovering novel neurobiological features that may be important for the recovery of acuity and new treatment targets.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.