0% READ
FULL TEXT
INTRODUCTION
Brain-Derived Neurotrophic Factor (BDNF) is a homodimeric member of the neurotrophin family that plays a fundamental role in neuronal differentiation, survival, synaptic plasticity, and long-term potentiation. Through activation of the Tropomyosin Receptor Kinase B (TrkB) receptor, BDNF orchestrates critical signaling cascades governing neurogenesis and synaptic remodeling (Colucci-D’Amato et al., 2020). Owing to its profound neuroprotective and neuroregenerative properties, recombinant BDNF has been investigated as a therapeutic candidate for a wide spectrum of neurological and psychiatric disorders, including Alzheimer’s disease, Parkinson’s disease, Huntington’s disease, spinal cord injury, major depressive disorder, and ischemic brain injury.
Despite its compelling therapeutic potential, the development of BDNF as a protein-based biologic is severely constrained by intrinsic aggregation liability (Majid and Khan, 2023). Aggregation of protein therapeutics is a major pharmaceutical challenge, as it can compromise conformational integrity, reduce receptor-binding affinity, diminish biological potency, accelerate chemical degradation, and increase immunogenic risk. From a regulatory perspective, aggregated species are considered critical quality attributes due to their potential to trigger adverse immune responses (Pham and Meng, 2020). Therefore, identifying and mechanistically characterizing aggregation-prone determinants within BDNF is essential for enabling rational stabilization strategies.
Protein aggregation is a complex, multi-factorial phenomenon governed by structural, energetic, and physicochemical parameters (Rojekar et al., 2025). Key drivers include exposure of hydrophobic residues, aromatic π-π stacking interactions, β-sheet nucleation propensity, conformational flexibility of glycine-rich segments, electrostatic imbalance, and local energetic destabilization (Rabe et al., 2009). Importantly, aggregation propensity is rarely uniform across a protein structure; rather, it is typically dominated by discrete spatially clustered hotspots that function as nucleation centers.
While numerous computational tools exist to predict aggregation-prone regions, most studies rely on a single algorithm or sequence-based predictors alone (Navarro and Ventura, 2022). However, aggregation is inherently multidimensional, and robust mechanistic inference requires integration of structure-based scoring, intrinsic sequence propensity, solvent accessibility analysis, and energetic decomposition. A statistically supported, chain-resolved structural interrogation is therefore necessary to define aggregation-driving nuclei with high predictive confidence.
In this study, we perform an integrated computational dissection of aggregation determinants in BDNF at residue-level resolution. By combining structure-informed aggregation scoring, sequence-driven aggregation prediction, β-sheet aggregation modeling, solvent accessibility quantification, energetic contribution analysis, and statistical distribution profiling across 293 residue entries, we delineate dominant aggregation nuclei and establish a mechanistic hierarchy of aggregation drivers. This integrative framework provides actionable insights for rational protein engineering and formulation optimization of BDNF as a therapeutic biologic.
MATERIALS AND METHODS
Structural Dataset and Residue-Level Data Extraction
The biologically active homodimeric structure of human Brain-Derived Neurotrophic Factor (BDNF) (PDB ID: 1BND; resolution 2.30 Å) was retrieved from the Protein Data Bank. Structural preparation was performed using the Protein Preparation Wizard in BioLuminate (Schrödinger Release 2023-4, Schrödinger, LLC, New York, NY, USA). Preparation steps included bond order assignment, addition of hydrogen atoms, correction of formal charges, optimization of hydrogen-bonding networks, and restrained minimization using the OPLS4 force field with heavy-atom convergence criteria of 0.3 Å RMSD.
A total of 293 residue-level entries across Chain A and Chain B were analyzed. For each residue, quantitative descriptors including AggScore-based aggregation score, Solvent-Accessible Surface Area (SASA), relative side-chain solvent accessibility ratio, patch surface area, and energetic contributions (positive destabilizing, negative stabilizing, and hydrophobic components) were extracted using the BioLuminate residue interaction analysis module.
Multi-Algorithm Aggregation Prediction
Aggregation propensity was evaluated using three complementary computational frameworks to ensure multidimensional assessment of aggregation determinants. AggScore analysis was performed within the BioLuminate environment using default structural aggregation parameters. Aggregation-prone residues were defined as those exceeding the 75th percentile of the residue-level aggregation score distribution. Surface hydrophobic clustering thresholds and spatial exposure criteria were retained at software-defined default settings to preserve methodological consistency.
Intrinsic sequence-based aggregation propensity was assessed using the online Aggrescan server (version 2023). Predictions were conducted using default parameters, including a window size of five residues and the experimentally derived a3v intrinsic aggregation propensity scale. Residues exhibiting positive normalized aggregation propensity values were classified as aggregation-prone segments.
To further evaluate β-sheet-driven aggregation potential, Zyggregator analysis was performed using its standard physicochemical weighting framework, which integrates hydrophobicity, charge distribution, and secondary structure propensity. Default algorithm coefficients were applied without modification. Regions with aggregation probability scores exceeding the server-defined threshold were considered aggregation-susceptible.
The outputs from AggScore, Aggrescan, and Zyggregator were subsequently integrated by cross-referencing residue-level predictions across all three methods to identify concordant aggregation hotspots. This multi-algorithm approach enabled simultaneous evaluation of structure-dependent exposure, intrinsic sequence-encoded aggregation propensity, and β-sheet-mediated aggregation mechanisms.
Energetic Contributions and Solvent Accessibility Analysis
Residue-wise energetic decomposition was performed using the BioLuminate module (Schrödinger, LLC, New York, NY, USA) following structural preparation of the BDNF dimer (PDB ID: 1BND). The prepared structure was subjected to residue-level interaction energy analysis to quantify individual energetic contributions within the protein microenvironment. For each residue, total interaction energy was decomposed into stabilizing (negative) and destabilizing (positive) components based on intramolecular interaction profiles. Additionally, hydrophobic energetic contributions were calculated in kcal/mol to estimate the extent of non-polar interactions contributing to local aggregation-driving forces. These energetic descriptors were extracted for all residues across both chains and incorporated into the statistical dataset.
Side-chain solvent accessibility was quantified using Solvent-Accessible Surface Area (SASA) calculations implemented within BioLuminate. Relative solvent accessibility ratios were computed by normalizing residue exposure against standard reference values for fully exposed amino acids (Tien et al., 2013). This normalization enabled identification of residues with elevated surface exposure. Because aggregation requires intermolecular contact, solvent-exposed hydrophobic and aromatic residues were considered particularly relevant for aggregation susceptibility.
To evaluate spatial clustering of aggregation determinants, aggregation-prone residues were mapped onto the protein surface, and patch surface areas were calculated based on contiguous solvent-accessible regions exceeding aggregation score thresholds. Patch Area (Ų) was determined by summing the surface contributions of neighbouring residues within defined spatial proximity (Jones and Thornton, 1997). This approach enabled identification of extended aggregation-competent domains rather than isolated single-residue vulnerabilities.
All energetic, accessibility, and patch parameters were exported as residue-level quantitative descriptors and subjected to subsequent statistical analysis and chain-wise comparisons.
Statistical Analysis
All residue-level quantitative descriptors were statistically evaluated to characterize aggregation behavior within the BDNF dimer. Descriptive statistics including mean, Standard Deviation (SD), minimum, maximum, and percentile distributions were calculated for aggregation scores, hydrophobic energetic contributions, stabilizing and destabilizing interaction energies, solvent accessibility ratios, and patch surface areas. To assess inter-variable relationships, Pearson correlation analysis was performed between aggregation score and hydrophobic contribution, aggregation score and solvent accessibility, and hydrophobic contribution and solvent accessibility. Correlation strength was interpreted based on Pearson’s correlation Coefficient (r), and linear regression analysis was conducted to quantify predictive relationships. Coefficient of Determination (R²) values were calculated to assess the proportion of variance explained. Comparative analysis between Chain A and Chain B was conducted using unpaired two-tailed t-tests to determine statistical differences in aggregation score, hydrophobic contribution, and solvent accessibility. Statistical significance was defined at p<0.05. All analyses were performed using GraphPad Prism (version 10.0; GraphPad Software, USA). Results are expressed as mean±SD unless otherwise specified.
RESULTS
Aggregation Distribution and Statistical Characterization
The aggregation landscape of BDNF exhibits pronounced heterogeneity and statistical asymmetry (Table 1). Across 293 residues, the mean aggregation score was 2.78±3.26, with a maximum value of 17.13. The aggregation score distribution demonstrated strong positive skewness (skewness coefficient=2.20), and Shapiro-Wilk testing confirmed significant deviation from normality (W=0.681, p<0.0001). The marked disparity between the mean and maximum values quantitatively supports a nucleation-dominated aggregation pattern characterized by rare high-intensity hotspots rather than uniform structural instability.
| Parameter | Mean±SD | Minimum | Maximum |
|---|---|---|---|
| Aggregation Score | 2.78±3.26 | 0.28 | 17.13 |
| Patch Surface Area (Ų) | 147.55±195.52 | 16 | 1039.8 |
| Hydrophobic Energy (kcal/mol) | 0.43±0.3 | 0 | 1.53 |
| Positive Energy Contribution | 360.44±369.61 | 8.14 | 1268.31 |
| Negative Energy Contribution | −0.65±0.77 | −2.28 | 1.78 |
| SASA (%) | 29.39±19.25 | 1 | 60 |
Patch-level analysis further reinforced this heterogeneity. Sixty aggregation-associated surface patches were identified, exhibiting a mean surface area of 147.55±195.52 Ų, with values ranging from 16.00 to 1039.80 Ų. The patch size distribution demonstrated pronounced positive skewness (skewness coefficient=2.34), indicating the presence of rare but disproportionately large aggregation-competent domains. Shapiro-Wilk testing confirmed non-normal distribution of patch sizes (W=0.64, p<0.0001). Several extended patches exceeded 700 Ų, suggesting that aggregation susceptibility is driven by clustered surface domains rather than isolated residue-level vulnerabilities.
Hydrophobic energetic contributions averaged 0.43±0.3 kcal/mol, with moderate dispersion across residues, whereas stabilizing (negative) energetic components demonstrated broader variability. This energetic heterogeneity suggests that local destabilization and hydrophobic clustering are unevenly distributed across the BDNF structure, reinforcing the concept that aggregation propensity is concentrated within discrete structural microenvironments.
Chain-wise Comparison of Aggregation and Energetic Parameters
Chain-wise analysis revealed differential aggregation behavior between the two monomeric subunits of BDNF. Residue-level aggregation scoring identified a dominant high-intensity hotspot in Chain A centered at GLY62, which exhibited the highest aggregation score (17.13) among all analyzed residues (Figure 1). This hotspot formed a contiguous cluster with PRO60, MET61, and TYR63, spanning residues 60-63. Collectively, these residues demonstrated elevated solvent accessibility and increased hydrophobic energetic contribution relative to neighboring regions. In addition, this segment displayed higher positive energetic contributions, indicative of localized destabilization consistent with aggregation-prone microdomains.
The elevated aggregation score at GLY62, together with the hydrophobic contribution of MET61 and the aromatic character of TYR63, suggests that this segment represents a structurally vulnerable region within Chain A (Figure 2). The spatial clustering of aggregation-prone residues within this narrow sequence window supports the presence of a localized aggregation nucleus rather than diffuse structural instability. Statistical validation further confirmed the exceptional nature of this region: GLY62 exhibited a Z-score of 4.71 relative to the global aggregation score distribution, confirming its classification as an extreme statistical outlier. The mean aggregation score for Chain A was 2.94±4.08, whereas Chain B exhibited a mean of 2.01±2.57. Independent two-tailed t-test analysis revealed a statistically significant difference between the two chains (t=2.08, p=0.039), indicating modest but statistically significant differences in aggregation score distribution between the two chains.
In contrast, Chain B exhibited a broader but comparatively lower-intensity aggregation profile. Moderate aggregation scores were observed for aromatic residues such as PHE56 and TYR55, consistent with their potential involvement in hydrophobic and π-π interactions. Additionally, the C-terminal segment spanning residues 95-110 displayed aggregation signatures across multiple predictive frameworks, including β-sheet propensity models, suggesting a secondary aggregation-prone region that may contribute to aggregation propagation rather than primary nucleation.
Correlation between Energetic Contributions and Solvent Accessibility
Residues identified as high-intensity aggregation hotspots exhibited elevated hydrophobic energetic contributions alongside reduced stabilizing (negative) energy components when compared with the overall residue population. These residues also demonstrated increased solvent accessibility ratios, indicating enhanced surface exposure (Figure 3). Pearson correlation analysis demonstrated a statistically significant but very weak negative association between hydrophobic energetic contribution and aggregation score (r=-0.129, p=0.039; R²=0.017). Although localized hotspot regions exhibited elevated hydrophobic contributions, the global residue-level dataset did not demonstrate a strong linear relationship, indicating that hydrophobicity alone is insufficient to predict aggregation propensity. This finding suggests that aggregation susceptibility is not driven solely by hydrophobic magnitude but is influenced by the combined effects of structural exposure, local energetic destabilization, and spatial clustering.
Notably, convergence of elevated aggregation score, increased solvent accessibility, and enhanced hydrophobic energetic contribution was observed within residue ranges 55-65 and 95-110. These segments consistently displayed higher aggregation-related parameters relative to surrounding regions, supporting their identification as dominant aggregation-prone domains within the BDNF structure.
DISCUSSION
The present integrated computational analysis indicates that aggregation susceptibility in BDNF is not uniformly distributed across the protein structure but is concentrated within specific residue clusters. Descriptive statistical analysis demonstrated a right-skewed distribution of aggregation scores, with a limited number of residues exhibiting markedly higher values than the overall mean. This observation suggests that aggregation propensity may be driven by localized high-intensity segments rather than generalized structural instability.
Chain-wise comparison revealed that Chain A contains a dominant aggregation hotspot within residues 60-63, particularly centered on GLY62. This region displayed the highest aggregation score along with increased hydrophobic energetic contribution and elevated solvent accessibility. The concurrence of these parameters suggests that this microdomain represents a structurally vulnerable segment prone to intermolecular association. The elevated hydrophobic contribution observed for MET61, together with the aromatic nature of TYR63 and the conformational characteristics of GLY62 and PRO60, may collectively enhance local aggregation susceptibility through surface exposure and hydrophobic interactions. Although glycine is intrinsically non-hydrophobic and does not directly promote hydrophobic aggregation, its conformational flexibility may facilitate transient structural fluctuations that expose adjacent hydrophobic residues such as MET61 and TYR63. Therefore, the aggregation propensity associated with GLY62 likely arises from its structural positioning within a flexible loop region rather than from intrinsic physicochemical aggregation potential. Structural mapping onto the BDNF crystal structure (PDB ID: 1BND) indicates that residues 60-63 are located within a loop region connecting β-strands. Loop regions are often associated with increased solvent exposure and structural mobility, which may enhance aggregation susceptibility under destabilizing conditions.
In contrast, Chain B exhibited a broader but lower-intensity aggregation profile. Although moderate aggregation scores were observed in aromatic residues such as PHE56 and TYR55, (Elhassan et al., 2025). The absence of an extreme high-magnitude hotspot comparable to GLY62 suggests quantitative differences in hotspot intensity between the two chains; however, given the homodimeric symmetry of BDNF, functional asymmetry cannot be concluded without experimental validation. The C-terminal segment (residues 95-110), which showed aggregation signatures across predictive models, may contribute to aggregation propagation rather than primary nucleation.
Given that BDNF exists as a homodimer, the spatial positioning of aggregation-prone residues relative to the dimer interface may influence aggregation behavior. Preliminary structural mapping suggests that the 60-63 segment appears spatially positioned in proximity to the dimer interface based on structural visualization, which may contribute to the quantitative differences in aggregation score distribution observed between Chain A and Chain B. The correlation between hydrophobic energetic contribution and aggregation score further emphasizes the importance of hydrophobic clustering in BDNF aggregation. Although localized hydrophobic enrichment was observed within hotspot regions, global statistical analysis indicates that hydrophobicity alone does not account for aggregation susceptibility. Residues with elevated hydrophobic energy values and higher solvent accessibility were consistently associated with increased aggregation propensity (Malleshappa Gowder et al., 2014). These findings indicate that surface-exposed hydrophobic regions are likely to facilitate intermolecular contacts, promoting aggregation under destabilizing conditions.
From a pharmaceutical perspective, identification of discrete aggregation-prone segments provides opportunities for targeted stabilization strategies (Rahban et al., 2023). Site-directed modification of residues within the 60-63 region could potentially reduce aggregation susceptibility; however, structural and functional validation would be necessary to ensure preservation of biological activity. In addition, formulation strategies incorporating hydrotropic agents or surfactants capable of masking hydrophobic surfaces may further enhance stability (Zoeller et al., 2022). Stabilization of the C-terminal aggregation-prone segment may also contribute to minimizing aggregation during storage and handling.
Overall, the integration of aggregation scoring, energetic contribution analysis, solvent accessibility assessment, and chain-wise comparison provides a coherent framework for understanding aggregation behavior in BDNF and supports rational approaches for improving its stability as a therapeutic protein. It is important to note that the present findings are derived from computational predictive models and residue-level energetic decomposition. While statistically robust, these results require experimental validation through biophysical aggregation assays, mutational analysis, or structural stability measurements to confirm their functional relevance.
CONCLUSION
The present multi-algorithm computational analysis suggests that aggregation propensity in BDNF is concentrated within specific residue clusters rather than being uniformly distributed across the protein structure. Chain-wise evaluation identified a prominent aggregation-prone segment within residues 60-63 of Chain A, characterized by high aggregation score, elevated hydrophobic energetic contribution, and increased solvent accessibility. Secondary aggregation-prone regions were observed within the C-terminal segment, which may contribute to aggregation progression.
The combined analysis of aggregation scoring, sequence-based prediction, energetic contribution assessment, and solvent accessibility profiling provides a comprehensive understanding of structural regions associated with aggregation susceptibility. These findings may support the rational design of stabilization strategies, including targeted residue modification and formulation approaches aimed at reducing aggregation and improving the stability of BDNF as a therapeutic protein. Future experimental studies will be necessary to validate the mechanistic implications of these computational findings and to confirm the functional relevance of the identified aggregation-prone regions.