Tunable Energetic Stabilization Governs Self- and Mixed-Aggregation of Möbius and Bracelet Cyclotides ()
1. Introduction
Conventional synthetic pesticides continue to dominate the global pest control market owing to their broad-spectrum activity, extensive distribution networks, and relatively low cost. Despite their effectiveness, the environmental and ecological consequences of these agents remain significant and increasingly concerning. For example, Glyphosate, the most widely used herbicide in the United States and the active ingredient in Roundup, has been associated with unintended effects on non-target organisms [1] [2]. These include chelation of essential mineral nutrients in plants and disruption of beneficial gut microbiota such as Snodgrassella alvi in bumble bees. Glyphosate exerts its herbicidal action by inhibiting 5-enolpyruvylshikimate-3-phosphate synthase (EPSPS), a key enzyme in the shikimate pathway responsible for the biosynthesis of the aromatic amino acids’ phenylalanine, tyrosine, and tryptophan [3]. Because this pathway is conserved among plants, bacteria, and fungi, glyphosate-based formulations may produce broad off-target effects beyond their intended use. These ecological and public health concerns have accelerated the search for safer, more sustainable pest control alternatives.
Among emerging bioactive compounds, plant-derived cyclotides have gained increasing attention as promising eco-friendly candidates. Interest in cyclotides originated in the 1960s following ethnobotanical observations of a traditional uterotonic tea prepared from Oldenlandia affinis in the Democratic Republic of the Congo [4] [5]. The active peptide responsible for this effect, later named Kalata B1, displayed remarkable stability even after prolonged boiling, leading to the identification of the first member of the cyclotide family. Cyclotides are small, backbone-cyclized peptides typically comprising 28 - 37 amino acid residues stabilized by a cystine knot motif that confers exceptional thermal, chemical, and proteolytic resistance [6]. They are primarily found in plants belonging to the Rubiaceae (coffee) and Violaceae (violet) families and are classified into two principal subfamilies—Möbius and Bracelet—based on differences in backbone topology and structural features. The Bracelet subfamily, characterized by the absence of a backbone twist, constitutes the majority of known cyclotides and exhibits greater structural diversity, whereas Kalata B1 belongs to the Möbius class [6].
Beyond their extraordinary stability, cyclotides display diverse biological activities, including anthelmintic, nematocidal, cytotoxic, insecticidal, antifouling, and antimicrobial effects [5]. These activities largely arise from their ability to disrupt cellular membranes through pore formation and membrane permeabilization, initiated by insertion of hydrophobic surface patches into lipid bilayers or membrane-mimetic micelles [7]. Additionally, Kalata B1 has been reported to enter mammalian cells through both endocytosis and direct membrane translocation [8]. Such properties have stimulated interest in their application as next-generation antimicrobials, where they have demonstrated significant activity against pathogens such as Acinetobacter baumannii. Collectively, these characteristics position cyclotides as attractive scaffolds for the development of sustainable biopesticides.
Despite extensive characterization of their biological functions, comparatively little is understood about the structural and physicochemical determinants governing cyclotide self-association and aggregation. Elucidating these interactions is important for understanding stability, bioavailability, and formulation behavior, all of which directly influence their practical deployment in agricultural settings. In this study, we investigate cyclotide self-association mechanisms using computational molecular dynamics simulations. We model interactions among identical and heterogeneous cyclotides across varying concentrations and analyze aggregate formation alongside potential, van der Waals, and electrostatic (Coulombic) energy contributions to identify the dominant forces driving assembly. By clarifying the molecular basis of cyclotide aggregation, this work aims to inform the rational design and optimization of cyclotide-based biopesticides for sustainable pest management.
2. Methodology
2.1. Cyclotide Structure Preparation
Cyclotide structures were obtained from the RCSB Protein Data Bank using their corresponding four-character PDB identifiers [9]. Five representative cyclotides were selected to capture sequence and structural diversity: Cycloviolacin O2 (CO2; 2KNM), Kalata B1 (1JJZ), Varv F (3E4H), Pase A (7K7X), and Circulin A (1BH4). Kalata B1, Circulin A, and Pase A belong to the Möbius subfamily, whereas Varv F and CO2 belong to the Bracelet subfamily (Figure 1). These structures were used to evaluate how differences in sequence composition and topology influence aggregation behavior.
2.2. Molecular Dynamics Simulations
All molecular dynamics (MD) simulations were performed using the ChemCompute GPU-accelerated platform, which utilizes Nanoscale Molecular Dynamics (NAMD) as the underlying simulation engine [10]. Simulations were executed on the Expanse GPU cluster configuration (10 processor cores, 91 GB RAM) with a maximum allowable runtime of 48 hours to ensure sufficient equilibration and convergence of aggregation systems [11].
Two aggregation regimes were modeled:
Self-aggregation systems, consisting exclusively of Kalata B1 (kB1) at varying peptide counts (2, 5, 10, and 20 molecules).
Mixed-aggregation systems, composed of Möbius and Bracelet cyclotides in equivalent total peptide counts.
Figure 1. Ribbon representations of the initial cyclotide structures used for molecular dynamics simulations. Shown are the three-dimensional biological assemblies arranged from left to right and top to bottom: (a) CO2, (b) kB1, (c) Varv F, (d) Pase A, and (e) Circulin. All structures were retrieved from the RCSB Protein Data Bank and prepared prior to simulation.
All systems were solvated in explicit TIP3P water within a cubic simulation box of 175 × 175 × 175 Å3. This box dimension was selected to minimize periodic boundary artifacts and prevent premature intermolecular interactions with periodic images, particularly in higher-concentration simulations. Simulations were conducted under constant volume and temperature (NVT ensemble) at 310 K. This temperature approximates warm agricultural environments and the physiological temperature range of common insect pests. Temperature regulation was achieved using Langevin dynamics to ensure stable thermal control throughout each trajectory. Each system was propagated for 1,000,000 integration steps using a 2 fs timestep. This timestep balances computational efficiency with numerical stability. The extended trajectory length was empirically selected to ensure that potential energy profiles reached stable plateaus indicative of aggregation equilibration.
2.3. Force Field and Treatment of Non-Bonded Interactions
Simulations employed established biomolecular force fields (CHARMM/AMBER families as implemented in NAMD), which describe the total potential energy of the system as the sum of bonded and non-bonded contributions. The AMBER99SB force field was used for these simulations. The total potential energy is expressed as:
(1)
Bonded interactions include bond stretching, angle bending, and torsional rotations. Non-bonded interactions consist of van der Waals (vdW) and electrostatic (Coulombic) terms (Equation (1)). van der Waals interactions were computed for all non-bonded atom pairs within a defined cutoff distance to ensure computational efficiency while preserving short-range dispersion interactions critical for peptide packing.
Long-range electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method. PME is an efficient implementation of Ewald summation that partitions electrostatic calculations into real-space and reciprocal-space components under periodic boundary conditions. This approach enables accurate computation of long-range Coulombic interactions without truncation artifacts.
By approximating the infinite summation of electrostatic interactions across periodic images, PME ensures reliable modeling of charge-charge interactions in solvated peptide systems. The use of PME was essential for accurately capturing screened electrostatics within explicit water environments and for maintaining energetic consistency across varying system sizes and compositions.
2.4. Trajectory and Statistical Analyses
Following completion of the simulations, time-resolved energy terms were extracted directly from the ChemCompute/NAMD trajectory log files. The following energetic quantities were analyzed:
Total potential energy (PE);
van der Waals (vdW) energy;
Coulombic (electrostatic) energy.
Total potential energy was used as a quantitative descriptor of aggregate thermodynamic stabilization. The vdW and Coulombic components were evaluated independently to determine the relative contributions of dispersion-driven packing interactions and screened electrostatic interactions in both self-aggregation and mixed-aggregation systems.
Energy-time datasets were exported and analyzed using Microsoft Excel. To assess equilibration and stability convergence, energies were compared across trajectories and aggregation conditions.
Statistical significance was evaluated using a two-factor analysis of variance (ANOVA). The two factors were:
1) Aggregation type (self-aggregation vs mixed-aggregation);
2) Aggregate concentration (2, 5, 10, and 20 starting cyclotide molecules).
A significant threshold of p < 0.05 was applied for all energy components (PE, vdW, and Coulombic).
3. Results
3.1. Aggregation Morphology from Molecular Dynamics Simulations
Molecular dynamics simulations were performed to compare the self-association of Kalata B1 (kB1) with mixed-composition cyclotide systems across increasing peptide numbers (2, 5, 10, and 20 molecules). Representative final configurations revealed concentration-dependent aggregation behavior in both cohorts, with distinct structural trends between self- and mixed-aggregation systems (Figure 2 and Figure 3). In the kB1 self-association simulations, increasing peptide number produced progressively more compact assemblies (Figure 3). The 2- and 5-peptide systems formed small clusters, whereas the 10- and 20-peptide systems exhibited dense, closely packed aggregates with minimal inter-peptide separation. This indicates enhanced intermolecular association at higher concentrations. In contrast, mixed cyclotide systems displayed fewer compact structures as concentration increased (Figure 2). Although lower-count systems formed cohesive clusters, larger systems adopted more dispersed or branched arrangements. Thus, aggregation compactness increased monotonically with concentration in the self-association cohort (Figure 3) but decreased in the mixed cohort, demonstrating that peptide composition influences aggregate morphology (Figure 2).
![]()
Figure 2. Representative aggregation configurations of mixed cyclotide systems obtained during molecular dynamics simulations. (a) Co-aggregation of kB1 and CO2. (b) Mixed aggregation comprising two kB1, two CO2, and one Circulin A molecules. (c) Aggregation of two kB1, two CO2, two Circulin A, two Pase A, and two Varv F molecules. (d) Larger-scale aggregation consisting of four kB1, four CO2, four Circulin A, four Pase A, and four Varv F molecules.
Figure 3. Representative aggregation states of kB1 cyclotides observed during molecular dynamics simulations at increasing system sizes: (a) two kB1 molecules, (b) five kB1 molecules, (c) ten kB1 molecules, and (d) twenty kB1 molecules.
3.2. Potential Energy Profiles and Thermodynamic Stabilization
To quantify aggregate stability, total potential energy was monitored throughout each trajectory. All systems exhibited a rapid initial decrease in potential energy followed by convergence to a stable plateau, indicating equilibration and attainment of thermodynamically stable configurations. Despite substantial differences in peptide composition and concentration, the final potential energies were highly similar across all simulations. Self-aggregation systems converged to approximately −1.54 to −1.58 × 106 kcal·mol−1, while mixed systems reached comparable endpoints within the same range. Statistical comparison of initial and final energies revealed significant differences between cohorts (p < 0.05). The similarity in converged energies indicates that both self- and mixed-aggregation processes result in assemblies of comparable thermodynamic stability (Figure 4); however, the subtle differences are impactful. Furthermore, equilibration occurred within the early portion of the trajectories, suggesting that stable intermolecular contacts form rapidly under the simulated conditions. Notably, systems exhibiting distinct morphologies (e.g., dispersed versus compact aggregates) nevertheless displayed nearly identical potential energies, indicating that structural compactness does not directly correlate with global energetic stability.
3.3. van der Waals Interaction Energies
van der Waals (vdW) interaction energies were evaluated to assess the contribution of short-range non-bonded contacts to aggregation. Initial vdW energies varied with system size and composition; however, all simulations converged to similar final values. Self-aggregation trajectories converged to approximately 2.25 - 2.33 × 105 kcal·mol−1, while mixed systems reached comparable values within the same interval (Figure 5). Convergence occurred rapidly, consistent with early establishment of favorable intermolecular packing interactions. These results indicate that dispersion-dominated contacts contribute similarly to stabilization across both homogeneous and heterogeneous cyclotide assemblies, despite differences in sequence composition.
Figure 4. Time evolution of potential energy during aggregation of systems composed of 2, 5, 10, and 20 molecules for (a) kB1 cyclotides and (b) mixed cyclotide systems.
Figure 5. Time evolution of van der Waals interaction energies during aggregation of systems composed of 2, 5, 10, and 20 molecules for (a) kB1 cyclotides and (b) mixed cyclotide systems.
3.4. Coulombic Interaction Energies
Electrostatic contributions were assessed through Coulombic interaction energies. As observed for vdW and total potential energies, initial values differed among systems but converged to nearly identical endpoints following equilibration. Final Coulombic energies ranged from approximately −2.09 × 106 to −2.14 × 106 kcal·mol−1 for both cohorts (Figure 6), with statistically significant differences observed between the groups (p < 0.05. The timescale of convergence closely matched that of the vdW and total energy profiles, indicating coordinated stabilization of electrostatic and short-range interactions during aggregate formation. The statistical significance of the differences in convergence energy across systems suggests that differences in residue charge distribution may slightly alter the overall energetic landscape of aggregation under the simulated conditions.
Across all simulations, three consistent features emerged. First, systems rapidly equilibrated, reaching stable energy plateaus within short simulation timescales. Second, the potential, van der Waals, and Coulombic energy components converged to comparable magnitudes yet statistically different, regardless of peptide composition, indicating similar underlying interaction strengths across systems. Third, distinct aggregate morphologies led to small but statistically significant differences in total energy. Collectively, these results suggest that cyclotide aggregation is governed by conserved intermolecular interaction networks that produce comparable thermodynamic stability across diverse compositions. Thus, peptide identity modulates aggregate geometry and packing arrangements and may play a role in association within the systems examined.
![]()
Figure 6. Time evolution of Coulombic interaction energies during aggregation of systems composed of 2, 5, 10, and 20 molecules for (a) kB1 cyclotides and (b) mixed cyclotide systems.
4. Discussions
The present molecular dynamics study provides a quantitative characterization of cyclotide self-association and mixed-composition aggregation at atomistic resolution. Across all simulated systems, aggregation was accompanied by rapid energetic stabilization and convergence of the potential, van der Waals, and Coulombic interaction energies to nearly identical magnitudes, despite substantial differences in peptide composition and concentration (Figures 4-6). These findings indicate that cyclotide association is governed primarily by conserved physicochemical interaction principles rather than sequence-specific effects.
4.1. Conserved Energetic Landscape of Cyclotide Aggregation
A central observation of this work is the convergence of total potential energies toward stable equilibrium plateaus for both homogeneous and heterogeneous assemblies. In molecular simulations, convergence to steady minima typically reflects occupation of well-defined basins on the free-energy surface. Although all systems reached thermodynamically stable states, statistical analysis revealed significant differences in the final equilibrium energies across aggregation types and concentrations (p < 0.05).
These results indicate that while cyclotides adopt broadly similar packing regimes, the corresponding free-energy basins are not energetically identical. Instead, sequence-dependent structural features introduce measurable differences in stabilization magnitude. This finding is notable given the structural diversity between members of the Möbius and Bracelet subfamilies. Variations in loop regions, residue composition, and charge distribution appear to modulate intermolecular interaction strengths rather than being completely averaged out at the aggregate scale.
The observed energetic differences likely arise from differential contributions of hydrophobic packing, dispersion interactions, and screened electrostatics in aqueous solution [12]. Thus, aggregation appears to be governed by shared physicochemical driving forces, but with stabilization amplitudes that depend on molecular composition and aggregate concentration.
The rapid equilibration observed across trajectories further supports this interpretation. Energy plateaus were reached during the early stages of simulation, indicating that favorable intermolecular contacts form efficiently and that large energetic barriers to association are absent. Such behavior is consistent with diffusion-limited aggregation driven primarily by short-range attractive forces, albeit modulated by sequence-specific interaction strengths.
4.2. Role of Non-Bonded Interactions
Decomposition of the energetic contributions provides additional mechanistic insight. Both van der Waals (vdW) and Coulombic energy components exhibited comparable convergence toward stable equilibrium plateaus across all systems. However, statistical analysis revealed significant differences in their final equilibrium magnitudes as a function of aggregation type and concentration (p < 0.05).
The general conservation of convergence behavior in vdW energies indicates that aggregate stabilization is predominantly driven by dispersion-mediated packing interactions. While the magnitude of vdW stabilization varied significantly between systems, the overall energetic trends suggest that dense intermolecular packing remains the principal thermodynamic driver. This is consistent with the compact aggregate structures observed in self-aggregation simulations and with the established tendency of cyclic peptides to minimize solvent-exposed hydrophobic surface area [12] [13].
Similarly, Coulombic energies approached stable equilibria across trajectories, though statistically significant differences were detected among systems. In aqueous environments, dielectric screening reduces the effective range of charge-charge interactions, limiting but not eliminating sequence-dependent electrostatic contributions. Thus, electrostatics appear to modulate aggregate stability rather than dictate fundamentally distinct aggregation pathways.
Taken together, these findings support a model in which cyclotide aggregation is governed primarily by non-specific short-range intermolecular contacts, with sequence-dependent variations influencing the magnitude—but not the fundamental mechanism—of stabilization [12] [14] [15].
4.3. Morphology-Energy Decoupling
Although overall energetic trends were conserved, distinct differences in aggregate morphology were observed between self- and mixed-composition systems. Self-association of Kalata B1 produced increasingly compact clusters with increasing concentration, whereas mixed systems formed more dispersed or branched assemblies at higher peptide numbers (Figure 2 and Figure 3).
While statistical analysis revealed significant differences in total equilibrium energies between systems (p < 0.05), these quantitative variations did not directly correlate with the observed morphological differences. In other words, aggregates exhibiting distinct geometries were found to occupy energetically similar regimes, despite measurable differences in stabilization magnitude.
This partial decoupling between morphology and thermodynamic stabilization suggests that multiple structural arrangements may reside within closely spaced regions of the free-energy landscape. Cyclotide assemblies therefore likely populate a range of metastable conformations that differ geometrically yet remain comparably favorable energetically [15]. Such degeneracy is common in soft-matter and peptide aggregation processes, where numerous packing configurations satisfy similar energetic constraints. Consequently, aggregate compactness alone may not serve as a reliable proxy for thermodynamic stability.
The energetic conservation observed here has practical implications for the design of cyclotide-based biotechnological or agricultural agents. If aggregation thermodynamics are largely sequence-independent, formulation strategies may not require strict control over peptide composition to achieve stable assemblies. Instead, properties such as bioactivity, target specificity, or environmental persistence can be optimized without substantially altering aggregation stability. Moreover, the rapid self-association observed in simulations suggests that cyclotides can readily form stable aggregates under physiologically relevant conditions, potentially enhancing resistance to degradation and prolonging functional lifetime. Such behavior may be advantageous for sustained delivery applications, including biopesticide or therapeutic formulations [16].
5. Conclusions
The broadly conserved aggregation behavior observed across all systems provides useful insight for the rational design of cyclotide-based biotechnological and agricultural agents. Although statistically significant differences in equilibrium energies were detected as a function of aggregation type and concentration (p < 0.05), all systems converged toward stable thermodynamic plateaus governed by similar intermolecular interaction profiles. This indicates that cyclotide aggregation is driven by shared physicochemical principles—primarily dispersion-mediated packing and screened electrostatics—while allowing sequence-dependent modulation of stabilization magnitude.
Such behavior suggests that stable assemblies can be achieved across diverse cyclotide compositions, providing flexibility for optimizing functional properties such as bioactivity, target specificity, or environmental persistence without fundamentally disrupting aggregation mechanisms. Rather than strict sequence-independence, the present results support a model of mechanistic conservation with quantitative energetic tunability. The rapid self-association observed in all simulations further indicates that cyclotides readily form compact and energetically favorable aggregates under physiologically relevant aqueous conditions. This intrinsic propensity for association may enhance resistance to degradation and extend functional lifetime, potentially benefiting sustained-delivery applications, including biopesticide and therapeutic formulations.
Despite these promising implications, several limitations of the present computational framework must be acknowledged. The simulations were conducted under idealized conditions that do not fully capture the environmental heterogeneity of real biological or agricultural systems, where fluctuations in temperature, pH, ionic strength, and interactions with membranes or soil constituents may significantly influence aggregation kinetics and stability. Additionally, molecular dynamics results remain dependent on the selected force-field parametrization, which approximates non-bonded interactions and may not fully account for long-range or many-body effects. Finally, the finite system sizes and simulation timescales accessible in this study may preclude observation of slower or larger-scale aggregation pathways.
Future investigations should therefore extend these simulations to larger ensembles and longer trajectories, incorporate heterogeneous environmental models, and integrate experimental validation to determine how the computationally identified energetic trends translate to practical applications.
Acknowledgements
This work was supported by the National Science Foundation (Award No. 2107567). The authors gratefully acknowledge Georgia Gwinnett College and its School of Science and Technology for financial and institutional support. The University of Yaoundé I, Department of Inorganic Chemistry, is also acknowledged for academic support. During the preparation of this manuscript, ChatGPT was used to assist with language editing and organizational improvements. All AI-assisted content was carefully reviewed and edited by the authors, who assume full responsibility for the accuracy, integrity, and originality of the final work.