Introduction

Polyketide Synthases (PKS) and Non-Ribosomal Peptide Synthetases (NRPS) are modular biological assembly lines responsible for many clinically important therapeutics. Traditional synthetic biology has sought to generate novel macrocycles by manually swapping mega-enzymes between pathways. However, as the MIBiG 4.0 database has expanded to over 46,000 constituent proteins, random combinatorial searches have become statistically intractable without systematic scoring heuristics.

We hypothesized that lightweight optimization methods—specifically sequential decision processes and thermodynamic optimization—could autonomously prune this massive sequence space and optimize structural junctions without large-scale neural network training. While deep-learning approaches offer powerful predictive capabilities, they introduce substantial dependency management overhead and reproducibility challenges that limit accessibility in resource-constrained settings.

Existing computational systems biology platforms generally fall into two categories. Database-centric exploration suites (e.g., ClusterCAD) provide excellent manual exploration interfaces but lack autonomous generative design logic. Conversely, deep-learning suites like DeepBGC and antiSMASH provide state-of-the-art predictive classification but require complex dependency stacks (PyTorch, TensorFlow) that can be difficult to deploy deterministically. GenerativeBGCs occupies a middle ground: a functional generative pipeline relying on core computing paradigms (sequential optimizers, thermodynamic schedules, alignment-free sequence comparison) within standard-library Python, with optional external validation endpoints for 3D structural assessment.

Results

Statistical Validation of AI-Guided Assembly

Our core finding, derived from a three-way ablation comparing Monte Carlo, UCB1 exploration, and full AI (UCB1 + SA) pipelines: the integration of structured optimization significantly improves interface compatibility compared to random assembly baselines.

1,000 bootstrap resamples on the generated Top-50 distribution yield clear advantages:

• Baseline (Monte Carlo): Mean SIS 96.95 [95% CI: 96.53–97.38]
• Full AI (MDP-inspired + SA): Mean SIS 97.54 [95% CI: 97.21–97.87]

A two-sided paired permutation test (10,000 permutations) confirms the performance delta is not attributable to random variation (Δ = +0.589, p < 0.0001). Cohen's d effect size is reported alongside p-values to characterize the practical magnitude of improvement. We note that this absolute improvement on a compressed SIS scale is modest; however, the consistency of the effect across bootstrap resamples and the statistical significance of the paired test support the conclusion that the optimization regimen provides a reliable, non-trivial advantage in boundary compatibility selection.

Hyperparameter Sensitivity & Distributional Robustness

A frequent concern in computational biology is high sensitivity to hyperparameter values. We swept the thermodynamic cooling constants ($T \in [0.8, 0.9]$) and UCB1 exploration parameters using actual chimera generation runs, confirming intrinsic robustness. Across all parameter configurations, the algorithmic outputs deviated by a bounded margin. This suggests that the optimization approach provides stable performance across a reasonable hyperparameter range, reducing the risk of overfitting to specific parameter choices.

Sequence Plausibility Assessment and Optional 3D Validation

Evaluating the top 10 GenerativeBGC candidates via our local k-mer Markov transition filter confirmed robust biological plausibility at the sequence level. As an optional supplementary validation step (not part of the core zero-dependency pipeline), selected candidate junction sequences can be submitted to the Meta ESMFold API via standard urllib requests. The retrieved coordinate structures (.pdb) provide independent structural plausibility assessment. We emphasize that this external validation step is orthogonal to the core generative logic and requires network access to the ESM Atlas service.

Discussion

This work demonstrates that structured optimization heuristics can provide meaningful improvements in combinatorial biosynthetic pathway design without requiring heavy compute clusters or complex deep-learning dependency stacks. By employing fundamental computing paradigms—bandit-based exploration, computational thermodynamics, and alignment-free sequence comparison—executed within Python's standard library, we provide a deterministic and reproducible generative pipeline.

We acknowledge several limitations. First, the composite SIS metric, while incorporating multiple biophysical properties (hydropathy, charge, aromatic packing, secondary structure propensity), remains a heuristic simplification; physical expression may face 3D steric hindrance, allosteric effects, and docking domain specificities not captured without atomistic molecular dynamics simulations or experimentally determined co-crystal structures. Second, while the k-mer Jaccard measure provides a useful zero-dependency sequence homology proxy, dedicated HMM-based tools (HMMER) or pairwise alignment methods (BLAST) would offer stronger evolutionary evidence for tailoring gene substitution. Third, the sequential decision process uses an ε-greedy policy without learned value functions; incorporating proper value function approximation or Monte Carlo tree search could further improve design quality.

GenerativeBGCs produces ready-to-synthesize .gbk sequences representing computationally plausible, biologically grounded synthetic operons, pending host-specific codon optimization and experimental validation. The implementation, verified by both bootstrap confidence intervals, permutation tests, and effect size reporting, represents a step toward accessible, reproducible computational tools for secondary metabolite engineering.

Methods

Cryptographic Data Verification

The pipeline executes solely on the local MIBiG 4.0 dataset. To preclude silent data drift, the execution environment enforces an SHA-256 hash lock (4b196343ed...).

Sequential Decision Process for Module Assembly

Pathway construction is formulated as a sequential decision process inspired by the MDP framework, defined by the tuple $(S, A, T, R)$:

• State $s_k \in S$: The partial chimeric assembly line at step $k$, represented as an ordered list of selected proteins $[p_0, p_1, \ldots, p_{k-1}]$.
• Action $a_k \in A$: A binary choice—retain the host protein at position $k$, or substitute a donor protein sampled from the donor pool.
• Transition $T(s_k, a_k) \rightarrow s_{k+1}$: Deterministic—the selected protein is appended to the assembly, advancing to state $s_{k+1}$.
• Reward $R(s_k, a_k)$: The composite Structural Interface Score (SIS) at the boundary between protein $p_{k-1}$ and the newly selected protein at position $k$.
• Policy $\pi$: An $\varepsilon$-greedy policy ($\varepsilon = 0.15$) that selects the action maximizing boundary SIS with probability $1 - \varepsilon$, and explores the alternative with probability $\varepsilon$.

We note that this formulation uses a fixed, non-learned policy rather than iterative policy improvement. The ε-greedy mechanism provides exploration-exploitation balance without requiring training episodes or value function approximation. Future work could incorporate learned value functions to improve decision quality.

Composite Structural Interface Score (SIS)

The SIS evaluates inter-protein boundary compatibility via four weighted biophysical components:

1. Hydropathy differential (weight 0.35): Absolute difference in mean Kyte-Doolittle hydropathy (Kyte & Doolittle 1982) between the C-terminal docking window (30 aa) of the upstream protein and the N-terminal window of the downstream protein.
2. Electrostatic charge compatibility (weight 0.25): Net charge product at physiological pH (~7.4), penalizing like-charge repulsion at the boundary interface.
3. Aromatic packing complementarity (weight 0.20): Differential in aromatic residue fraction (F, W, Y, H), reflecting interface packing compatibility.
4. Secondary structure propensity mismatch (weight 0.20): Chou-Fasman helix propensity differential (Chou & Fasman 1978) combined with Gly/Pro flexibility index asymmetry.

The composite penalty is converted to a compatibility score via Gaussian decay: $\text{SIS} = 100 \cdot e^{-\gamma \cdot \text{penalty}^2}$, where $\gamma$ is auto-calibrated from natural intra-cluster boundaries to yield a mean SIS of ~90 for known-functional interfaces.

Thermodynamic Simulated Annealing (SA)

Poorly compatible structural boundaries (SIS < 70) are candidates for computational linker insertion. Rather than greedy selection, the platform employs Simulated Annealing: suboptimal conformations are accepted with decaying probability governed by the Boltzmann distribution ($e^{-\Delta E / T}$), reducing the risk of premature convergence to local optima. The cooling schedule uses $T_0 = 5.0$ with multiplicative decay rate 0.85.

Two-Tier Tailoring Gene Substitution

Secondary metabolites rely on downstream auxiliary genes (e.g., methyltransferases, oxidoreductases). Gene substitution employs a two-tier matching strategy:

1. Primary: K-mer Jaccard sequence similarity — An alignment-free measure computing the Jaccard coefficient between overlapping tri-peptide (k=3) sets of the host and donor auxiliary protein sequences. Candidates exceeding a Jaccard threshold of 0.25 are considered homologous. We acknowledge that while k-mer Jaccard provides a useful lightweight proxy, dedicated profile HMM methods (e.g., HMMER) would provide more sensitive detection of remote homologs.
2. Secondary: TF-IDF annotation similarity — When multiple donor candidates have similar k-mer scores (within 0.05 of the maximum), functional annotation text is compared using TF-IDF cosine similarity as a tiebreaker (tokenization dropping fragments ≤ 2 characters, minimum cosine threshold 0.40). This leverages existing MIBiG metadata to disambiguate among sequence-similar candidates.

Offline K-Mer Markov Evaluation and Optional ESMFold Validation

GenerativeBGCs implements a di-peptide Markov chain evaluator within the standard library to verify base sequence plausibility. As an optional external validation step (not required for the core pipeline), junction sequences can be submitted to the Meta ESMFold API via zero-dependency urllib requests. This provides supplementary 3D structural plausibility assessment but introduces a dependency on an external service. We clearly distinguish this optional validation from the self-contained core generative pipeline.

Data and code availability

Pipeline code (GenerativeBGC main generation and orchestration engine) and the statistical ablation suite: https://github.com/yzjie6/GenerativeBGCs. Data: MIBiG 4.0 sequence database. External validation tools: antiSMASH API, ESMFold API (optional). The core generative pipeline is fully self-contained within Python's standard library.

References

• Auer, P., Cesa-Bianchi, N., & Fischer, P. (2002). Finite-time analysis of the multiarmed bandit problem. Machine learning, 47(2), 235-256.
• Blin, K., Shaw, S., Kloosterman, A. M., Charlop-Powers, Z., van Wezel, G. P., Medema, M. H., & Weber, T. (2019). antiSMASH 5.0: updates to the secondary metabolite genome mining pipeline. Nucleic acids research, 47(W1), W81-W87.
• Chou, P. Y., & Fasman, G. D. (1978). Prediction of the secondary structure of proteins from their amino acid sequence. Advances in enzymology and related areas of molecular biology, 47, 45-148.
• Hannigan, G. D., Prihoda, D., Palicka, A., Soukup, J., Klempir, O., Rampula, L., ... & Medema, M. H. (2019). A deep learning genome-mining strategy for biosynthetic gene cluster prediction. Nucleic acids research, 47(18), e110-e110.
• Kirkpatrick, S., Gelatt, C. D., & Vecchi, M. P. (1983). Optimization by simulated annealing. Science, 220(4598), 671-680.
• Kyte, J., & Doolittle, R. F. (1982). A simple method for displaying the hydropathic character of a protein. Journal of molecular biology, 157(1), 105-132.
• Medema, M. H., Kottmann, R., Yilmaz, P., Cummings, M., Biggins, J. B., Blin, K., ... & Glöckner, F. O. (2015). Minimum information about a biosynthetic gene cluster. Nature chemical biology, 11(9), 625-631.
• Terlouw, B. R., Blin, K., Navarro-Muñoz, J. C., Avalon, N. E., Chevrette, M. G., Egbert, S., ... & Medema, M. H. (2023). MIBiG 3.0: a community-driven effort to annotate biosynthetically active genomic regions. Nucleic acids research, 51(D1), D603-D610.