\documentclass[11pt]{article}
\usepackage{amsmath,amsthm,amssymb,amscd}
\usepackage{hyperref}
\usepackage{geometry}
\usepackage{tikz-cd}
\usepackage{booktabs} % For better table aesthetics
\usepackage{graphicx}
\usepackage{subcaption} % provides the subfigure environment

\usepackage{multirow}
\usepackage{tikz}
\usetikzlibrary{calc,positioning,3d}
\usepackage[normalem]{ulem} % cross through
%\usepackage{makeidx}
%\makeindex
\usepackage{imakeidx}
\usepackage{makecell}
\usepackage{float}
\makeindex
\usepackage{seqsplit}
\newcommand{\ttseq}[1]{{\ttfamily\seqsplit{#1}}}
\usepackage{xr-hyper}
\externaldocument{SI}    
\usepackage{cleveref}   

%\usepackage{subfigure} % For subfigures

\geometry{margin=1in}
\parindent=0pt
\parskip=6pt

\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{remark}[theorem]{Remark}


\title{GBNL: Graded Betti Number Learning of Complex Biological Data}

\author{Mushal Zia$^1$,  
	Faisal Suwayyid $^{1,2}$, 
	and Guo-Wei Wei\footnote{Corresponding author: Guo-Wei Wei (weig@msu.edu).~}$^{~1,3,4}$ \\		
	$^1$Department of Mathematics,\\
	Michigan State University, MI 48824, USA.\\
	$^2$Department of Mathematics,\\
	King Fahd University of Petroleum and Minerals, Dhahran 31261, KSA.\\
	$^3$Department of Mathematics and Statistics,\\
	Michigan State University, MI 48824, USA.\\
	$^4$Department of Biochemistry and Molecular Biology,\\
	Michigan State University, MI 48824, USA.
}


\date{} 
\begin{document}
	\maketitle 
	
	\begin{abstract}	    
While persistent homology is widely used for data shape analysis, persistent commutative algebra (PCA) has seen limited adoption in machine learning and data science. Unlike persistent homology, which delivers topological invariants in the form of Betti numbers, PCA provides both algebraic invariants and graded Betti numbers. However, graded Betti numbers have seldom been applied to real-world data. In this work, we introduce the first-of-its-kind application of commutative algebra graded Betti numbers in machine learning and data science. Specifically, we present Graded Betti Number Learning (GBNL) for protein-nucleic acid binding prediction. Protein-DNA/RNA interactions are fundamental to cellular processes such as replication, transcription, translation, and gene regulation, and their understanding and prediction remain challenging. GBNL represents each nucleic acid sequence as a family of $k$-mer-specific sets and derives persistent graded Betti invariants from PCA, generating multiscale topological representations of local nucleotide organization. To incorporate cross-molecule context, these graded Betti representations are paired with transformer-based protein embeddings, linking nucleotide-level signals with global protein patterns. The proposed graded Betti representations effectively detect single-site mutations and distinguish complete mutation patterns. Operating on primary sequences with minimal preprocessing, GBNL bridges commutative algebra, reduced algebraic topology, combinatorics, and machine learning, establishing a new paradigm for comparative sequence analysis. Numerical studies using three datasets highlight the success of GBNL in protein-nucleic acid binding prediction.
		
	\end{abstract}
	
	Keywords: Persistent commutative algebra, graded Betti numbers, Stanley-Reisner theory, machine learning, protein-nucleic acid binding.   
	
	\newpage	
	{\setcounter{tocdepth}{4} \tableofcontents}
	\setcounter{page}{1}
	\newpage	
	
	\section{Introduction}
	
	Protein-nucleic acid interactions underpin the molecular framework of gene regulation and flow of genetic information. Through sequence-specific recognition of DNA and RNA, proteins orchestrate replication, repair, transcription, splicing, and translation, thereby coupling genotype to controlled cellular function. Moreover, the interpretation and execution of genetic instructions inside a cell are correctly governed by the specificity and stability of these contacts. Quantification of protein-nucleic acid binding affinities is thus critical for not only understanding recognition principles but also for enabling synthetic biology, genome engineering, and therapeutic innovation~\cite{bertoldo2023rna}. At atomic resolution, affinity emerges from a calibrated balance of hydrogen bonds, electrostatics, hydrophobic packing, and shape complementarity. Disruptions in any of these processes can miswire regulation and are linked to neurodegeneration, autoimmunity, and cancer~\cite{kisby2024introducing}. Therefore, conserving both global organization and nucleotide-resolution details remains a persistent challenge molecular biophysics and computational biology.
	
	Standard binding assays such as electrophoretic mobility shift assays (EMSA), surface plasmon resonance (SPR), isothermal titration calorimetry (ITC), fluorescence readouts, and filter binding, deliver accurate affinities yet they remain costly and scale poorly~\cite{stockley2009filter}. Computational substitutes reduce cost however, they introduce tradeoffs. Moreover, physics-based free energy calculations such as thermodynamic integration, free energy perturbation, and molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) yield detailed estimates but become prohibitive for large or flexible complexes. Furthermore, knowledge-based potentials \cite{zhang2005knowledge}, empirical and force-field scoring functions % 	\cite{nithin2019structure,yin2008medusascore}
	\cite{nithin2019structure}, and machine learning with engineered descriptors \cite{bitencourt2019machine} improve speed but can lack transferability and interpretability. These approaches perform well for protein-protein and protein-ligand systems, yet often degrade on DNA due to nucleic acid conformational variability and limited curated affinities. For RNA, crosslinking immunoprecipitation sequencing \cite{hafner2010transcriptome}, sequence driven predictors \cite{zhao2011structure,liu2016prediction,shen2023svsbi}, knowledge-based scoring \cite{tuszynska2011dars}, and coarse grained docking \cite{setny2011coarse} broaden coverage, but quantitative affinity measurements remain sparse because of molecular flexibility and experimental constraints \cite{iwakiri2016improved}.
	
	Recent progress in computational molecular biology has been driven by data and machine learning models that now play a central role in predicting biomolecular interactions and accelerating drug discovery \cite{lin2023evolutionary,song2024multiobjective}. The integration of bioinformatics with modern learning frameworks has enabled large scale inference of complex physical and chemical relationships from high-dimensional molecular data \cite{lo2018machine}. Within this evolution, mathematical artificial intelligence (AI), particularly topological deep learning (TDL) introduced in 2017 \cite{cang2017topologynet}, has emerged as a rational approach for representing and learning biomolecular structure \cite{papamarkou2024position,nguyen2020review} and has demonstrated strong performance in community challenges such as the D3R Grand Challenge series \cite{nguyen2019mathematical,nguyen2020mathdl}. This framework has since inspired the development of topological sequence analysis for genomic modeling \cite{hozumi2024revealing}, where delta complex formulations allow efficient treatment of large sequence data \cite{liu2025topological2}, and category theoretic formulations improve resolution among closely related genetic variants \cite{liu2025topological}. In parallel, machine learning models for protein-DNA binding have evolved from early atomic pairwise statistical potentials \cite{zhao2010structure}, sequence-based affinity profiling \cite{rastogi2018accurate}, and physics-informed methods such as DNAffinity \cite{barissi2022dnaffinity} to ensemble architectures including PreDBA \cite{yang2020predba}, PDA-Pred \cite{harini2023pda}, and emPDBA \cite{yang2023empdba}. Following the first quantitative protein-RNA affinity dataset~\cite{yang2013dataset}, the trajectory moved toward structure-guided prediction with methods by Deng et al.~\cite{deng2019predprba} and Nithin et al.~\cite{nithin2019structure}, with subsequent refinements by Hong et al.~\cite{hong2023updated} and Harini et al.~\cite{harini2024pred}. However, despite these advances, most approaches continue to rely on stacked regressors, class-specific architectures, or interface-level descriptors to mitigate data scarcity, underscoring the need for frameworks that combine mathematical rigor with biological interpretability.
	
	Persistent homology (PH) serves as a foundational method in TDL for quantifying multiscale shape analysis via the birth and death of topological features. The resulting representations built from Betti numbers, barcodes, and Betti curves yield concise encodings of protein structure, molecular surfaces, and genomic organization~\cite{nguyen2020review,wee2025review}. Because it records only topological invariants, PH can miss local geometric or physicochemical changes. The persistent Laplacian (PL), a spectral approach introduced in 2019 \cite{wang2020persistent,chen2019evolutionary}, augments PH with non-harmonic spectra to incorporate geometric details and has progressed in theory and applications, including protein engineering and viral variant analysis \cite{memoli2022persistent,qiu2023persistent,chen2022persistent}. Nevertheless, both PH and spectral approaches primarily emphasize global information and connectivity.
	
	Commutative algebra offers a complementary route by translating simplicial and combinatorial data into polynomial ideals and algebraic invariants \cite{eisenbud2013commutative}. Persistent Stanley-Reisner theory integrates multiscale analysis with this algebraic setting to produce persistent graded Betti numbers, persistent facet ideals, $f$-vectors, $h$-vectors, as algebraic analogues of topological summaries \cite{suwayyid2025persistent}. These new representations follow the evolution of square-free monomial ideals under filtration, revealing how geometric and physicochemical interactions emerge or vanish across scales. The same algebraic machinery has already improved structure based protein-ligand affinity prediction through a commutative algebra representation of interfacial geometry and chemistry \cite{feng2025caml}. Extending this idea, a commutative algebra neural network learns graded algebraic signatures that trace mutation driven changes and relate them to disease mechanisms \cite{wee2025commutative}. In parallel, a $k$-mer formulation casts sequence positions into persistent algebraic features, enabling comparative genomics and phylogenetic inference at scales \cite{suwayyid2025cakl}. Together, these directions establish persistent commutative algebraic invariants as interpretable, data efficient descriptors that connect topology, geometry, and learning for biomolecular systems.
	
	In this study, we present GBNL, the first graded Betti framework integrating commutative algebra with biomolecular modeling and machine learning. It encodes nucleic acid organization through graded Betti numbers tracked over a positional filtration for protein-nucleic acid binding prediction. We pair these algebraic/topological descriptors with transformer-based protein embeddings to capture interface chemistry at two complementary scales. The result is a compact, interpretable feature set for predicting protein-DNA and protein-RNA binding. We validate the approach on three standard datasets and a targeted mutation analysis of the N-China-F primer of the SARS-CoV-2 N gene, showing that the graded Betti signal isolates local sequence changes that conventional Betti counts miss. 
	
	
\section{Results and Discussion}\label{sec:results}

\subsection{Overview of the Workflow}

\begin{figure}[t!]
	\centering
	\begin{subfigure}[b]{1\textwidth}
		\includegraphics[width=\textwidth]{workflow.pdf}
	\end{subfigure}
	\caption{\textbf{An Illustration of GBNL workflow.}
		(a). Representative protein-nucleic acid complex provide biological context; downstream computations use only the DNA/RNA sequence (for per-nucleotide point clouds) and the protein amino-acid sequence (for embeddings).
		(b). \(k\)-mers are extracted from the DNA sequences of the given complexes. For each \(k=1\), the set of its occurrence positions within the sequence is treated as input data.
		(c). For each nucleotide, a Vietoris-Rips (VR) complex with maximum dimension 1 is constructed from its position-based point cloud across an \(\varepsilon\)-filtration. At each selected \(\varepsilon\), the maximal simplices define a simplicial complex whose Stanley-Reisner (SR) ideal \(I^{\varepsilon}\) is formed. The minimal free resolution of \(S/I^{\varepsilon}\), where \(S\) is a polynomial ring, is then computed in Macaulay2 (M2), yielding the graded Betti numbers \(\beta_{i,j}\), which serve as the topological features for the corresponding nucleotide. 
		(d). From the protein amino acid sequences of the same complexes, 2560-dimensional embedding vectors are generated using the state-of-the-art ESM2 model with 36 layers.
		(e). The extracted graded Betti features of DNA and the ESM2 embeddings of proteins are concatenated to form a unified algebraic-sequence representation.
		(f). Finally, a gradient boosted decision tree (GBDT) model is trained using these combined features for the prediction of experimental binding affinities.}
	\label{fig:workflow}
\end{figure}

Figure~\ref{fig:workflow} illustrates the overall workflow of the proposed GBNL framework for the prediction of protein-DNA/RNA binding. This approach couples algebraic invariants extracted from nucleic acid sequences with transformer-based representations of protein sequences, then learns a mapping from this combined representation to experimental binding affinities.

For each nucleotide, we construct position-based point clouds over a filtration and build Vietoris-Rips (VR) complexes with maximum simplex dimension 1. At each level, the complex is mapped to its Stanley-Reisner ideal and a minimal graded free resolution is computed to read off the graded Betti numbers \(\beta_{i,j}\). These \(\beta_{i,j}\) keys act as persistent algebraic signatures of local sequence organization and provide a compact, interpretable summary of multiscale combinatorial structure at the nucleotide level. In parallel, protein sequences are embedded with a modern transformer to capture global sequence regularities that complement the nucleotide specific algebraic features. The two feature streams are concatenated and used by a gradient boosted decision tree learning model to produce affinity predictions.

This design provides two benefits. First, the multiscale commutative algebra graded Betti numbers localize signal to specific nucleotides over filtration scales, which supports interpretation of model outputs. Second, the joint representation improves predictive accuracy by utilizing the rich transformer-based protein embeddings. Further details of protein sequence embedding, nucleic acid sequence processing, the datasets used and model parametrization are provided in the Methods section. 

\subsection{Protein-DNA/RNA Acid Binding Affinity Prediction}

\begin{figure}[t!]
	\centering
	\begin{subfigure}[b]{1\textwidth}
		\includegraphics[width=\textwidth]{datasets.pdf}
	\end{subfigure}
	\caption{\textbf{Prediction summary and dataset statistics for S186, S142, and S322.}
		(a). Normalized experimental vs.\ predicted binding affinities; each dataset is min-max scaled to its own unit range and offset for visual separation, with marginal density plots shown along the top and right.
		(b). Sequence-length distribution: dots indicate the observed count at each length, and the curves show a 5-nt centered rolling average.
		(c). Experimental Gibbs free energy ($\Delta G$) vs.\ sequence length: dots mark the per-length mean $\Delta G$ values, and the curves show the corresponding 5-nt rolling average.}
	\label{fig:datasets}
\end{figure}

Protein-nucleic acid recognition drives catalysis, transport, signaling, transcription, and translation, while preserving chromosome stability and regulating gene expression. Disruption of these contacts is linked to autoimmune disease, inherited disorders, and cancer. Binding affinity reflects the combined effects of hydrogen bonding, van der Waals contacts, dipole interactions, electrostatics, and hydrophobic forces, and analyzing their joint impact enables rational, structure based therapeutic design. In this study, we use a graded Betti framework from commutative algebra to encode sequence derived determinants of these interactions and to predict protein-nucleic acid binding affinity with an interpretable, data driven model. Using this representation, we predict protein-DNA/RNA acid binding affinities with a model designed to balance accuracy with mechanistic insight.

We first assess predictive accuracy against the SVSBI reference~\cite{shen2023svsbi}. The SVSBI model integrates an ESM-based transformer for proteins with DNABERT for nucleic acids and reports an average Pearson of 0.669 and an RMSE of 1.98\,kcal/mol. On the other hand, our GBNL framework achieves \(R_p=0.691\) with RMSE \(=1.81\)\,kcal/mol on S186 using 3829 features, surpassing SVSBI on both correlation and error as shown in Table~\ref{table:Rp1}. Performance remains stable on S142 and S322, with 4287 and 4567 feature vectors, respectively, with \(R_p=0.66\), RMSE \(=2.17\)\,kcal/mol and \(R_p=0.657\), RMSE \(=2.03\)\,kcal/mol (Table~\ref{table:Rp2}). These results indicate consistent generalization across datasets of different sizes and compositions.

\autoref{fig:datasets}a compares experimental and predicted binding affinities after min-max normalization within each dataset, indicating stable calibration and no systematic bias between predictions and measurements. In \autoref{fig:datasets}b, sequence-length distributions are depicted with observed counts and smoothed 5-nt rolling averages. S186 exhibits a tight length profile centered near 12$\sim$nt with most sequences confined between 10-18$\sim$nt and dropping sharply beyond 20$\sim$nt. S142 is dominated by short fragments below 15$\sim$nt (median 9$\sim$nt) and showing fewer occurrences of longer strands. In contrast, S322 exhibits a broader mid-length representation, maintaining substantial counts across the 8-17$\sim$nt interval before gradually tapering at longer lengths, pointing to greater structural heterogeneity among complexes.


\begin{table}[htb!]
	\small
	\centering
	\caption{
		Comparison of prediction performance between existing SVSBI model and our GBNL framework on S186 for protein-nucleic acid binding affinity prediction. Reported metrics include Pearson correlation coefficient (\(R_p\)) and RMSE values in kcal/mol. All results are averaged over twenty independent runs with different random seeds, and the average metric values are reported.}
	\label{table:Rp1}
	\vspace{10pt}   
	\begin{tabular}{l | c | c}
		\hline
		\textbf{Model} & \textbf{\(R_p\)} & \textbf{RMSE (kcal/mol)} \\
		\hline
		\textbf{GBNL}   & \textbf{0.691} & \textbf{1.81} \\
		SVSBI \cite{shen2023svsbi}   & 0.669  &  1.98\footnotemark \\ 
		\hline
	\end{tabular}
\end{table}

\footnotetext{%
	The original RMSE of 1.45 reported by Shen et al. \cite{shen2023svsbi} was not converted into kcal/mol; here we apply the factor 1.3633 to obtain 1.98 kcal/mol.}

To illustrate how binding free energy (\(\Delta G\)) varies with sequence length across the datasets, \autoref{fig:datasets}c depicts the corresponding experimental and predicted trends. For S186, experimental energies span from \(-15.0\) to \(-5.8\)\,kcal/mol (mean \(-9.73 \pm 1.87\)\,kcal/mol), whereas predictions fall between \(-13.5\) and \(-5.8\)\,kcal/mol (mean \(-9.79 \pm 1.44\)\,kcal/mol). Both experimental and predicted profiles show a clear tendency toward more negative affinities as strand length increases, reaching minima near 42$\sim$nt around \(-12.6\)\,kcal/mol. A comparable pattern emerges in S142, where experimental \(\Delta G\) values range from \(-16.9\) to \(-4.3\)\,kcal/mol (mean \(-9.48 \pm 2.11\)\,kcal/mol), while predictions occupy a narrower interval of \(-13.3\) to \(-6.1\)\,kcal/mol (mean \(-9.48 \pm 1.46\)\,kcal/mol). Despite the strong short-length bias observed in its distribution, longer sequences exhibit more stable binding, with experimental and predicted minima near 80 and 93$\sim$nt at \(-15.9\) and \(-13.27\)\,kcal/mol, respectively. 

\begin{figure}[t!]
	\centering
	\begin{subfigure}[b]{1\textwidth}
		\includegraphics[width=\textwidth]{ACGU.pdf}
	\end{subfigure}
	\caption{\textbf{Persistent graded Betti curves evolution across nucleotides for the RNA sequence \texttt{ACAUGUUUUCUGUGAAAACGGAG}.}
		Four panels correspond to nucleotides A, C, G, and U, showing how graded Betti counts vary with the filtration threshold~\(\varepsilon\). Curve groups are organized by homological index \(i\) and internal degree \(j\), as indicated in the legends.}
	\label{fig:ACGU}
\end{figure}


Furthermore, S322 follows the same overall trend, with experimental energies extending from \(-16.9\) to \(-4.3\)\,kcal/mol (mean \(-9.64 \pm 1.99\)\,kcal/mol) and predicted values spanning \(-13.4\) to \(-5.7\)\,kcal/mol (mean \(-9.62 \pm 1.46\)\,kcal/mol). Here too, both curves display increasingly negative \(\Delta G\) values with sequence length, culminating near 80–87$\sim$nt at roughly \(-13\) to \(-15.9\)\,kcal/mol. Collectively, these patterns indicate that stronger (more negative) experimental binding affinities are consistently associated with longer strands, even though the majority of sequences cluster at shorter lengths. Overall, the graded Betti framework aligns with the observed sequence-affinity relationships and maintains consistent calibration across datasets of different size and composition. On S186 in particular, our GBNL method delivers improved correlation and reduced error relative to the SVSBI baseline. A comprehensive comparison of experimental and predicted binding affinities for S186, S142, and S322 are presented in \autoref{tab:S186}, \autoref{tab:S142}, and \autoref{tab:S322}, respectively, along with the corresponding PDBIDs.


\begin{table}[htb!]
	\small
	\centering
	\caption{Prediction performance of GBNL on S142 and S322 datasets.}
	\label{table:Rp2}
	\vspace{8pt}
	\begin{tabular}{l | c | c}
		\hline
		\textbf{Dataset} & \textbf{\(R_p\)}  & \textbf{RMSE (kcal/mol)} \\
		\hline
		S142  & 0.66 & 2.17 \\
		S322 & 0.657 & 2.03 \\
		\hline
	\end{tabular}
\end{table}

To complement the quantitative evaluation and illustrate how graded features evolve across nucleotides, we visualize in \autoref{fig:ACGU} the persistent graded Betti curves in the GBNL framework for the RNA sequence \ttseq{ACAUGUUUUCUGUGAAAACGGAG} from S142. Each panel corresponds to a nucleotide (A, C, G, U) and shows how graded Betti counts vary with the filtration parameter~\(\varepsilon\), grouped by homological degree \(i\) and internal degree \(j\). Early curve separations and plateau differences reveal distinct algebraic organization among nucleotides, particularly in higher \((j-i)\) gaps, where broader dynamic ranges correspond to locally denser configurations. This visualization highlights the interpretive granularity of the graded Betti representation beyond aggregate performance metrics.

	
	\section{Discussion}\label{sec:Discussion}
	
	\subsection{Performance Analysis}
Building on these results, we find that Macaulay2 (M2) is generally efficient for graded Betti computations and scales well on higher performance cloud computing (HPCC) systems. However, performance degrades on longer or denser sequences because the associated Stanley-Reisner ideals develop very large generating sets and high syzygy complexity. An illustrative case in S186 is the protein-DNA complex 4LJR, which contains a 35$\sim$nt mononucleotide run of T. Under our position-based embedding, we mark the positions of a chosen token (e.g., all T’s) as 1-based indices along the sequence, and for each \(\varepsilon\) draw an edge between two positions if their index distance \(\le \varepsilon\). At \(\varepsilon=34\) these T positions become fully connected, so the VR(1) 1-skeleton is the complete graph \(K_{35}\) (here \(K_n\) denotes the complete graph on \(n\) vertices). For a graph-only complex, every triangle of \(K_{35}\) is a minimal nonface, hence a cubic generator of the Stanley-Reisner ideal; there are \(\binom{35}{3}=6545\) such generators in a single state. Resolving \(k[\Delta^{\epsilon}]=S/I^{\varepsilon}\) with thousands of degree-3 generators across many plateau endpoints in \(\varepsilon\) consumes time and memory. Computing maximum simplex dimension 2 is even more demanding because 4 or more cliques that are not promoted to 3-simplices become degree-4 and higher generators, accelerating the growth in the resolution.

To handle such cases while preserving the graded Betti output for graphs, we replace the computer algebra step with an exact combinatorial procedure specialized to VR(1). For each \(\varepsilon\), we set \(K=\lfloor \varepsilon \rfloor\) and work with the threshold graph that connects two indices when their gap is at most \(K\). Filtration events occur only at integer gaps, so Betti values are constant on plateaus and we evaluate at their right endpoints. In a VR(dim\(=1\)) complex, Hochster’s formula implies that only two diagonals of the graded Betti table can be nonzero, i.e., the disconnection diagonal \(\beta_{i,i+1}(k[\Delta(\varepsilon)])\), which counts how many \(j\)-vertex choices split into more than one component, and the cycle diagonal \(\beta_{i,i+2}(k[\Delta(\varepsilon)])\), which counts loops across all \(j\)-vertex choices,  additionally \(\beta_{0,0}=1\). We compute these exactly by direct counting: a small dynamic program yields the induced-subgraph component distribution for \(\beta_{i,i+1}\), and a closed identity combines total edge counts at threshold \(K\) with that distribution to obtain \(\beta_{i,i+2}\). Thus, the combinatorial evaluation reproduces the graded Betti tables for VR(1) while running time in seconds. We adopt this route only when M2 becomes prohibitively slow or explodes; for all other sequences we use M2, which remains efficient in our pipeline and benefits from HPCC using multiple jobs and subjobs.

We now place our GBNL findings in the context of prior sequence based predictors of protein-nucleic acid binding affinity. Existing models for protein-DNA \cite{zhao2010structure,yang2020predba,harini2023pda,yang2023empdba} and protein-RNA \cite{deng2019predprba,hong2023updated,harini2024pred} prediction have expanded the landscape. However, many depend on extensive hand-crafted descriptors, task specific tuning, and varied datasets, prohibiting direct comparison. Recent protein-DNA predictors range from potentials to ensembles. For example, emPDBA \cite{yang2023empdba} reports \(R_p=0.66\) with \(\mathrm{MAE}=1.24\,\mathrm{kcal/mol}\) after DNA subclassing (vs. \(R_p=0.12\), \(\mathrm{MAE}=1.64\,\mathrm{kcal/mol}\) unclassified). PDA-Pred \cite{harini2023pda} achieves \(R_p=0.78\) with categorization (vs. \(R_p=0.21\) without class categorization). On a 36-complex blind set, emPDBA attains \(R_p=0.53\) with \(\mathrm{MAE}=1.11\,\mathrm{kcal/mol}\); PreDBA \cite{yang2020predba} \(R_p=0.30\); DDNA3 \cite{zhao2010structure} \(R_p=0.09\) with \(\mathrm{MAE}=1.80\,\mathrm{kcal/mol}\). 

On the other hand, protein-RNA models are fewer. PRA-Pred \cite{harini2024pred} trains on 217 complexes clustered at 25\% identity, partitions RNA into five structural groups and proteins by function, and with jack-knife selection achieves \(R_p=0.77\) with \(\mathrm{MAE}=1.02\,\mathrm{kcal/mol}\); on the standard 44-complex blind set it reports \(R_p=0.60\) with \(\mathrm{MAE}=1.47\,\mathrm{kcal/mol}\). Similarly, PRdeltaGPred \cite{hong2023updated} reaches \(R_p=0.41\) with \(\mathrm{MAE}=1.83\,\mathrm{kcal/mol}\) while PredPRBA \cite{deng2019predprba} yields \(R_p=0.07\) with \(\mathrm{MAE}=2.07\,\mathrm{kcal/mol}\) (best subclass \(R_p=0.48\)). This demonstrates that comparison across literature is complicated by heterogeneous sequence identity thresholds and different subclass protocols utilized in each study. Moreover, for many methods, there are no reported codes, datasets, or protocol, making it difficult to reproduce the claimed results.  

In contrast, our GBNL framework uses curated, biochemically consistent datasets, including S186, S142 and a cleanly annotated extension, S322. It encodes paired sequences via algebraic grades and indices across multiscale filtrations and learns with a single gradient boosting regressor. Operating directly on primary sequences, the method requires no three-dimensional structures, subclass labels, or handcrafted interface features, yielding a single unified predictor for protein-DNA and protein-RNA binding affinity. Instead of enumerating hydrogen bonds, base-pair numbers, or stacking motifs, the GBNL framework summarizes interactions at the sequence level, avoiding reliance on subclass-specific models and enabling high-throughput analysis independent of curated annotations.

\begin{figure}[h!]
	\centering
	\begin{subfigure}[b]{1\textwidth}
		\includegraphics[width=\textwidth]{covid_1.pdf}
	\end{subfigure}
	\caption{\textbf{Single mutation analysis by GBNL for N\text{-}China\text{-}F.}
		Rows: nucleotides. Column 1: PH Betti curves; columns 2 and 3: persistent graded Betti, all shown before and after mutation (MUT).}
	\label{fig:covid1}
\end{figure}

	\subsection{Model Interpretability}
	
Topological frameworks such as persistent homology (PH) is a well-established and powerful tool or multiscale quantification of biomolecular data. Yet, although PH captures global shifts in connectivity and cycles,it can obscure the underlying algebraic structures associated with specific mutations. Unlike PH summaries, GBNL separates homological content into graded components indexed by \(i\) and indices \(j\). Each \(\beta_{i,j}\) count generators (syzygies) at homological grade \(i\) and internal degree \(j\), yielding a precise account of how local perturbations shape the filtration landscape. Consequently, graded Betti curves uncovers intrinsic, nontrivial alterations that standard Betti summaries overlook.

We demonstrate in \autoref{fig:covid1} how GBNL sensitively captures a single-base mutation in the N-China-F primer of the SARS-CoV-2 N gene. A single nucleotide replacement in the reference sequence \ttseq{GGGGAACTTCTCCTGCTAGAAT} to \ttseq{AGGGAACTTCTCCTGCTAGAAT} at position 28881 (G\(\to\)A), one of the most frequent reported events (\(n=8672\)) \cite{wang2020mutations}, induces discernible shifts in topological and algebraic structure. At \(\varepsilon=0\), the A-specific component count increases from \(5\to6\) while loops remain zero. With larger \(\varepsilon\), PH shows modest rises in \(\beta_1\) (from \(2\to3\) at \(\varepsilon=13\) and \(6\to10\) at \(\varepsilon=20\)), indicating localized reorganization. Graded Betti profiles magnify these shifts: at the largest \(\varepsilon\), the gap \((j-i=2)\) for A increases \((1,3):10\to20\), \((2,4):15\to45\), \((3,5):6\to36\), while the corresponding G entries decrease \((1,3):20\to10\), \((2,4):45\to15\), \((3,5):36\to6\). Taken together, these opposing patterns reflect reallocation of generators and attenuation of higher-order syzygies near the mutated locus, aligning with a localized, single-site disturbance of the combinatorial neighborhood. PH and graded Betti curves for nucleotides C and T show no appreciable shift with mutation, confirming that all detected variation is attributable to the G\(\to\)A substitution. This analysis highlights the framework’s high sensitivity, showing that a single nucleotide alteration produces a measurable algebraic signature reflected in the graded Betti representation.

We next demonstrate how GBNL can efficiently capture the complete mutation in the same primer within the 22-base window. In \autoref{fig:covid_all}, the reference \ttseq{GGGGAACTTCTCCTGCTAGAAT} becomes \ttseq{AACGAATTTTTCTTGGTACAAT} \cite{wang2020mutations} after complete mutation. The triplet at 28881, 28882, and 28883 shows the highest counts (8672, 8660, 8659), and additional changes at 28887 (C\(\to\)T, 95), 28896 (C\(\to\)G, 52), 28899 (G\(\to\)C, 28), and 28900 (C\(\to\)T, 45) further perturb local connectivity. PH exhibits only limited variation in components and loops. Graded curves separate early in \(\varepsilon\), maintain broad plateaus, and span larger ranges at higher grades. For example, the component count shifts \(6\to9\) at \(\varepsilon=0\), and \(\beta_1\) increases from \(6\to20\) at \(\varepsilon=8\) and \(10\to28\) at \(\varepsilon=15\) for nucleotide T. At the largest \(\varepsilon\), the two-step gaps in the GBNL framework grow sharply: \((1,3):20\to84\), \((2,4):45\to378\), \((3,5):36\to756\), \((4,6):10\to840\). These values indicate a proliferation of algebraic generators across grades and scales. The graded Betti representation therefore reveals both precise single-site effects and the amplified structure induced by clustered substitutions. These cumulative changes illustrate how multiple concurrent mutations induce broad algebraic reorganization. The graded Betti representation exposing mutation-driven structural amplification far beyond what is visible through conventional persistent homology.

\begin{figure}[h!]
	\centering
	\begin{subfigure}[b]{1\textwidth}
		\includegraphics[width=\textwidth]{covid_all.pdf}
	\end{subfigure}
	\caption{\textbf{Complete mutation analysis by GBNL for N\text{-}China\text{-}F.}
		Rows: nucleotides. Column 1: PH Betti curves; columns 2 and 3: persistent graded Betti, all shown before and after mutation (MUT).}
	\label{fig:covid_all}
\end{figure}

Together, the single- and complete-mutation analyses highlight the dual interpretive strengths of the GBNL framework: it resolves fine-grained single-site changes and it captures the algebraic amplification that arises from multiple, co-occurring substitutions within the same genomic window. This interpretability has direct modeling value. By encoding geometric scale through \(\varepsilon\) and algebraic organization through homological degree \(i\) and internal index \(j\), the graded Betti feature space links mutation-driven combinatorial rearrangements to functional outcomes. The clear, grade-resolved separations between reference and mutated trajectories in the graded columns of \autoref{fig:covid1} and \autoref{fig:covid_all} align with the observed gains in predictive accuracy when graded features are used instead of descriptors based only on persistent homology. Consequently, graded Betti numbers serve as mathematically rigorous commutative algebraic invariants and biologically interpretable indicators of mutation-driven reorganization in biomolecular sequences.


	\section{Methods}
	\label{sec:methods}
	
	In this section, we describe the datasets, the vectorization pipeline based on persistent commutative algebra, and transformer-based featurization via ESM2 protein embeddings. Next, we present the persistent Stanley-Reisner framework and show how it yields the $k$-mer algebraic representation of sequences. Finally, we summarize the machine learning models and evaluation metrics. 
	
	\subsection{Datasets} 
	
	In this study, we use three datasets. Dataset~S186, introduced by Shen et al.~\cite{shen2023svsbi}, contains 186 protein-nucleic acid complexes with strand lengths from 5 to 45$\sim$nt. It was assembled from PDBbind-v2020 under a strict curation protocol: a complex was retained only if it had one unique protein sequence and one unique nucleic acid sequence; multiple chains sharing the same sequence were allowed, but any chain with mixed or hybrid bases (for example, both T and U) was discarded. We kept only entries with experimentally measured binding affinities at 298~K, removed records with ambiguous labels (\textasciitilde, $<$, or $>$), and excluded complexes with strands shorter than five nucleotides. The second dataset, S142, is a protein-RNA collection created by merging our curated PDBbind-v2020 entries with selected PRA-Pred~\cite{harini2024pred} complexes and keeping only RNA-binding proteins; RNA strand lengths range from 4 to 93$\sim$nt. Finally, S322 is a combined protein-nucleic acid set constructed in the same manner, with strand lengths from 5 to 93$\sim$nt. Together, these curated data collections span from short to long strands and support the analysis of sequence-dependent binding affinities.
	
	\subsection{Persistent Commutative Algebra Featurization for Nucleic Acids}\label{sec:Vectorization}
	
	For each nucleic acid sequence, we convert per-nucleotide information into a fixed-length vector of graded Betti features using a uniform, discrete filtration. The occurrence positions of the four mononucleotides (A, C, G, T/U) are treated as one-dimensional point clouds and evaluated on an integer \(\varepsilon\)-grid. At each \(\varepsilon\), we build a Vietoris-Rips (VR) complex on the corresponding positions for a maximum simplicial dimension 1, form the associated Stanley-Reisner ideal, and compute the minimal free resolution to obtain graded Betti numbers \(\beta_{i,j}\). This yields, for every nucleotide, a sequence of \(\beta_{i,j}\) values aligned across identical filtration scales, providing uniformly sampled algebraic summaries of sequence organization.
	
	Each graded Betti index \((i,j)\) corresponds to a distinct position in the Betti table, characterized by its homological degree \(i\) and internal degree \(j\). Consistent feature alignment is implemented by constructing the feature vector from the global union of all nonzero \((i,j)\) indices observed over the \(\varepsilon\)-grid. By default, we exclude \(\beta_{0,0}\). We populate features using a left-endpoint carry-forward scheme: at each \(\varepsilon\), the entry records the nearest preceding (or equal) nonzero \(\beta_{i,j}\) to preserve consistent feature evolution across the \(\varepsilon\)-grid. For example, if for some nucleotide~C, we observe \(\beta_{1,2}=3\) and \(\beta_{2,3}=2\) at \(\varepsilon=0\), then \(\beta_{1,2}=2\) and \(\beta_{2,3}=1\) at \(\varepsilon=7\), and finally \(\beta_{1,2}=1\) at \(\varepsilon=9\), the entries corresponding to \(\beta_{1,2}\) across \(\varepsilon=0,\dots,9\) become \([3,3,3,3,3,3,3,2,2,1]\), while those for \(\beta_{2,3}\) become \([2,2,2,2,2,2,2,1,1,0]\). Stacking the nucleotide-specific blocks produces a unified representation of graded Betti changes across the filtration scales with an identical schema for every sequence.
	
	For S186, we observe that the global unions comprise 25, 29, 27, and 67 distinct \((i,j)\) keys for A, C, G, and T/U, summing to 148. Given ten filtration steps (\(\varepsilon = 0\ldots9\)), each sequence obtains a nucleic acid feature vector of length 148 × 10 = 1480. In S142, the per-nucleotide key counts are \(\{A:49,\,C:51,\,G:61,\,U:31\}\), yielding a 1920-dimensional RNA feature vector. In S322, the counts are \(\{A:49,\,C:51,\,G:61,\,T/U:67\}\), yielding a 2280-dimensional nucleic acid feature vector. These representations provide concise algebraic summaries of sequence structure while preserving key organizational features across datasets.
	
	\subsection{Transformer-Based Protein Language Model}\label{sec:NLP}
	
	Recent advances in natural language processing have introduced transformer architectures capable of capturing contextual dependencies directly from primary protein sequences. We employ the Evolutionary Scale Modeling (ESM2) framework~\cite{lin2023evolutionary}, a state-of-the-art protein language model trained on large-scale sequence corpora to learn residue-level contextual embeddings without structural supervision. ESM2 encodes each amino acid sequence through 36 transformer layers and produces 2560-dimensional embeddings that summarize long-range biochemical and evolutionary interactions.
	
	These protein embeddings complement the persistent commutative algebra descriptors of nucleic acids by providing a parallel, sequence-driven characterization of the interacting protein partner. For each protein-nucleic acid complex, the ESM2 protein embedding is concatenated with the corresponding nucleic-acid feature vector to form the final model input. The resulting combined feature dimensions are 3829 for S186, 4287 for S142, and 4567 for S322 after removing features that are identically zero across all samples. These final design matrices are employed as the inputs for model training and evaluation. \autoref{fig:workflow} provides a schematic overview of the feature-generation pipeline.
	
	\subsection{Persistent Stanley-Reisner Theory}
	\label{subsec:PSRT}
	
	Our framework leverages persistent Stanley-Reisner theory (PSRT) to ground data analysis in commutative algebra. Whereas traditional TDA uses persistent homology to reveal geometric and topological patterns like loops and voids~\cite{su2025topological}, PSRT instead examines the underlying algebraic and combinatorial essence of simplicial complexes. Input data are initially mapped onto a simplicial complex simplicial complex of vertices, edges, triangles, and higher-order simplices to retain their topological and combinatorial structure. A filtration is then applied to monitor the emergence and persistence of these characteristics across multiple spatial or geometric scales, producing commutative-algebraic invariants such as persistent $h$-vectors, $f$-vectors, graded Betti numbers, and facet ideals~\cite{suwayyid2025persistent}. This study concentrates solely on graded Betti numbers, situating the investigation within an algebra-driven framework for data analysis.
	
	\subsubsection{Persistent Stanley-Reisner Structures Over a Filtration}
	
	Let \(k\) be a field, and let \(\Delta\) be a simplicial complex on the finite vertex set \(V = \{x_1, \dots, x_n\}\). Suppose \(f: \Delta \to \mathbb{R}\) is a monotone function, i.e., \(f(\tau) \le f(\sigma)\) whenever \(\tau \subseteq \sigma\), which induces an increasing filtration \(\Delta^{s} \subseteq \Delta^{\epsilon}\) for \(s \le \varepsilon\). Let \(S = k[x_1, \dots, x_n]\) be the standard graded polynomial ring over \(k\), and for each \(\varepsilon \in \mathbb{R}\), define the {Stanley-Reisner ideal} of \(\Delta^{\epsilon}\) as
	\[
	I^{\varepsilon} := \left\langle x_{i_1} \cdots x_{i_r} \,\middle|\, \{x_{i_1}, \dots, x_{i_r}\} \notin \Delta^{\epsilon} \right\rangle \subseteq S,
	\]
	with corresponding {Stanley-Reisner ring}
	\[
	k[\Delta^{\epsilon}] := S / I^{\varepsilon}.
	\]
	
	Since the filtration is increasing, the subcomplexes satisfy \(\Delta^{s} \subseteq \Delta^{\epsilon}\) for \(s \le \varepsilon\), which implies a descending chain of monomial ideals:
	\[
	I^s \supseteq I^{\varepsilon} \quad \text{for all } s \le \varepsilon.
	\]
	
	
	\subsubsection{Persistent Graded Betti Numbers of Stanley-Reisner Rings}
	
	Let \(k\) be a field and \(S = k[x_1, \dots, x_n]\) the standard graded polynomial ring. For each filtration level \(\varepsilon \in \mathbb{R}\), the Stanley--Reisner ring \(k[\Delta^{\epsilon}] := S / I^{\varepsilon}\) inherits a natural \(\mathbb{Z}\)-graded \(S\)-module structure and admits a minimal graded free resolution:
	\begin{equation}\label{eq:min-free-res}
		\cdots \longrightarrow
		\bigoplus_{j} S(-j)^{\beta_{i,j}(k[\Delta^{\epsilon}])}
		\longrightarrow \cdots \longrightarrow
		k[\Delta^{\epsilon}] \longrightarrow 0,
	\end{equation}
	where \(\beta_{i,j}(k[\Delta^{\epsilon}]) := \dim_k \operatorname{Tor}^S_i(k[\Delta^{\epsilon}], k)_j\) are the {graded Betti numbers}.
	
	Hochster’s formula relates these graded Betti numbers to the topological Betti numbers of the induced subcomplexes:
	\begin{equation}\label{eq:Hochster-general}
		\beta_{i,j+i}(k[\Delta^{\epsilon}])
		=
		\sum_{\substack{W \subseteq V \\ |W| = j+i}} 
		\dim_k \widetilde{H}_{j-1}(\Delta_W^{\varepsilon}; k),
	\end{equation}
	where \(\widetilde{H}_{j-1}(\Delta_W^{\varepsilon}; k)\) denotes the \((j-1)\)-st reduced simplicial homology group over \(k\), and 
	\(\Delta_W^{\varepsilon} := \{ \sigma \in \Delta^{\varepsilon} \mid \sigma \subseteq W \}\) is the subcomplex induced on the vertex set \(W \subseteq V\).
	
	In particular, Hochster’s formula can be reformulated in terms of the (non-reduced) Betti numbers of induced subcomplexes. For each integer \(i \ge 0\), the following identities hold:
	\begin{align}
		\beta_{i,i+1}(k[\Delta^{\varepsilon}]) &= \sum_{\substack{W \subseteq V \\ |W| = i+1}} \left( \beta_0(\Delta_W^{\varepsilon}) - 1 \right), \label{eq:Hochster-j1} \\
		\beta_{i,i+j}(k[\Delta^{\varepsilon}]) &= \sum_{\substack{W \subseteq V \\ |W| = i+j}} \beta_{j-1}(\Delta_W^{\varepsilon}), \quad \text{for all } j \ge 2, \label{eq:Hochster-j2}
	\end{align}
	where \(\Delta_W^{\varepsilon}\) denotes the subcomplex of \(\Delta^{\varepsilon}\) induced on the vertex subset \(W \subseteq V\), and \(\beta_r(\Delta_W^{\varepsilon})\) denotes the \(r\)-th Betti number of \(\Delta_W^{\varepsilon}\) with coefficients in \(k\). 
	
	To refine this in a persistent setting, for \(\varepsilon \le \varepsilon'\), we define the {persistent graded Betti number}
	\begin{equation}\label{eq:persistent-Betti-short}
		\beta_{i, i+j}^{\varepsilon,\varepsilon'}(k[\Delta]) 
		:= \sum_{\substack{W \subseteq V \\ |W| = i + j}}
		\dim_k \left( \iota_{j-1}^{\varepsilon,\varepsilon'} : \widetilde{H}_{j-1}(\Delta_W^\varepsilon) 
		\to 
		\widetilde{H}_{j-1}(\Delta_W^{\varepsilon'}) \right),
	\end{equation}
	where \(\iota_{j-1}^{\varepsilon,\varepsilon'}\) is the homomorphism on reduced homology induced by inclusion. This provides a multigraded algebraic refinement of classical persistent Betti numbers, encoding both topological persistence and the combinatorial properties of the evolving homology classes.
	
	In the special case where \(|W| = |V|\), the persistent graded Betti number reduces to
	\[
	\beta_{i, |V|}^{\varepsilon,\varepsilon'} = \beta_{|V| - i - 1}^{\varepsilon,\varepsilon'},
	\]
	recovering the classical persistent Betti number of homological degree \(|V| - i - 1\). More generally, the family \(\{\beta_{i, i+j}^{\varepsilon,\varepsilon'}\}_{i,j}\) encodes a richer multiscale invariant that interpolates between algebraic and topological persistence.
	
	
	Additional structural identities among the persistent Betti numbers further simplify this formula. In particular, it is known that
	\[
	\beta^{\varepsilon,\varepsilon'}_{0,0} = 1, \qquad
	\beta^{\varepsilon,\varepsilon'}_{i,i} = 0 \quad \text{for all } i \ge 1, \qquad
	\beta^{\varepsilon,\varepsilon'}_{0,j} = 0 \quad \text{for all } j \ge 1, \qquad
	\beta^{\varepsilon,\varepsilon'}_{i,j} = 0 \quad \text{for all } i > j.
	\]
	Consequently, one obtains
	\[
	B_0 = \beta^{\varepsilon,\varepsilon'}_{0,0} = 1,
	\]
	and for each \( j \ge 1 \), the alternating sum simplifies to
	\[
	B_j = \sum_{i=1}^{j-1} (-1)^i \beta^{\varepsilon,\varepsilon'}_{i,j}.
	\]
	
	
	\subsubsection{\(k\)-mer Algebraic Representations of Sequences}
	\label{subsec:kmer_representation}
	
	In this section, we specialize the \(k\)-mer algebra framework to the case \(k=1\). The framework, introduced by Hozumi et al.~\cite{hozumi2024revealing}, provides a principled approach for embedding sequences as collections of integer sequences in a geometric space. Through the Stanley-Reisner construction, this embedding induces an algebraic structure determined by the positional distribution of individual nucleotides.
	
	Let \(\mathcal{A}\) be a finite alphabet. A \(1\)-mer over \(\mathcal{A}\) is an element \(\boldsymbol{x}=x_1\in\mathcal{A}\).
	Given a fixed \(1\)-mer \(\boldsymbol{x}\), define the indicator function
	\[
	\delta_{\boldsymbol{x}}:\mathcal{A}\to\{0,1\},
	\qquad
	\delta_{\boldsymbol{x}}(\boldsymbol{y})=
	\begin{cases}
		1,& \boldsymbol{y}=\boldsymbol{x},\\
		0,& \boldsymbol{y}\neq\boldsymbol{x}.
	\end{cases}
	\]
	For a sequence \(S=s_1s_2\cdots s_N\in\mathcal{A}^N\), the set of starting positions at which \(\boldsymbol{x}\) occurs is
	\[
	S^{\boldsymbol{x}}
	=
	\Bigl\{\, i\in[1,N] \ \Bigm|\ 
	\delta_{\boldsymbol{x}}\bigl(s_i\bigr)=1 \Bigr\}.
	\]
	We regard \(S^{\boldsymbol{x}}\subset\mathbb{R}\) as a one–dimensional point cloud and define the pairwise distance matrix
	\[
	D^{\boldsymbol{x}}=\bigl(d^{\boldsymbol{x}}_{ij}\bigr)_{i,j\in S^{\boldsymbol{x}}},
	\qquad
	d^{\boldsymbol{x}}_{ij}=|i-j|.
	\]
	
	These distance matrices serve as inputs to compute algebraic features over a filtration interval \([r_0,r_1]\subset\mathbb{R}_{\ge 0}\).
	For \(r,r'\in[r_0,r_1]\) with \(r\le r'\), build the Vietoris--Rips complex on \(S^{\boldsymbol{x}}\) and denote by
	\[
	v^{r,r'}_{\boldsymbol{x}}=\bigl( v^{r,r'}_{i,j}(\boldsymbol{x}) \bigr)_{i<j\in\mathbb{Z}_{\geq 0}}
	\]
	the vector of algebraic invariants in each dimension \(i\) and internal degree \(j\).
	On the diagonal \(r=r'\), write
	\[
	v_{\boldsymbol{x}}=\bigl( v_{i,j}(\boldsymbol{x}) \bigr)_{i<j\in\mathbb{Z}_{\geq 0}}
	\]
	to represent the curves for $r\in [r_0,r_1]$.
	The full \(1\)-mer representation of \(S\) is obtained by concatenation over all \(\boldsymbol{x}\in\mathcal{A}\):
	\[
	\boldsymbol{v}^1_S
	:=
	\bigl( v_{\boldsymbol{x}} \ \bigm|\ \boldsymbol{x}\in\mathcal{A} \bigr)
	=
	\bigl( v_{i,j}(\boldsymbol{x}) \ \bigm|\ \boldsymbol{x}\in\mathcal{A} \bigr)_{i<j\in\mathbb{Z}_{\geq 0}},
	\]
	which we refer to as the \(1\)-mer algebraic representation of \(S\).
	
	In this study, we fix the nucleotide alphabet 
	\(\mathcal{A} := \{A,\, C,\, G,\, T/U\}\).
	Accordingly, the \(1\)-mer algebraic representation of a sequence \(S\) is given by
	\[
	\boldsymbol{v}^1_S
	=
	\bigl(
	v_{\boldsymbol{A}},\,
	v_{\boldsymbol{C}},\,
	v_{\boldsymbol{G}},\,
	v_{\boldsymbol{T/U}}
	\bigr),
	\]
	where each component \(v_{\boldsymbol{x}}\) encodes the algebraic invariants associated with the positional distribution of the nucleotide \(\boldsymbol{x}\in\mathcal{A}\).
	
	
	\subsection{Examples}
	
	\begin{figure}[t!]
		\centering
		% reduce spacing between subfigures by removing \hfill and adjusting width
		\begin{subfigure}[b]{0.24\textwidth}
			\includegraphics[width=\textwidth]{tetra.pdf}
		\end{subfigure}
		\hspace{6em} % small space between images
		\begin{subfigure}[b]{0.24\textwidth}
			\includegraphics[width=\textwidth]{bipyramid.pdf}
		\end{subfigure}
		
		\caption{\textbf{Geometric realizations of the example complexes.}
			(a). Tetrahedron with one missing face. (b). Square bipyramid (octahedron boundary) with each base vertex connected to both apices.}
		
		\label{fig:examples}
	\end{figure}
	
	\begin{example}[Tetrahedron with one missing face]\label{example1}
		Let $\Delta$ be the simplicial complex on $V=\{0,1,2,3\}$ obtained from a tetrahedron by removing the base face $\{1,2,3\}$. 
		All six edges are present,
		\[
		\{0,1\},\ \{0,2\},\ \{0,3\},\ \{1,2\},\ \{1,3\},\ \{2,3\},
		\]
		and the $2$-faces (facets) are exactly
		\[
		\{0,1,2\},\ \{0,1,3\},\ \{0,2,3\}.
		\]
		
		
		Let $S=k[x_0,x_1,x_2,x_3]$ and write $I_\Delta\subset S$ for the Stanley--Reisner ideal.
		The only minimal nonface is $\{1,2,3\}$, so $I_\Delta=(x_1x_2x_3)$ is a principal (squarefree) cubic, and the minimal free resolution of $S/I_\Delta$ is
		\[
		0 \longrightarrow 
		\underbrace{S(-3)}_{i=1}
		\longrightarrow 
		\underbrace{S}_{i=0}
		\longrightarrow S/I_\Delta \longrightarrow 0.
		\]
		By the standard dictionary (each copy of $S(-j)$ in homological degree $i$ contributes one to $\beta_{i,j}$), we obtain
		\[
		\beta_{0,0}=1,\qquad \beta_{1,3}=1,
		\]
		and $\beta_{i,j}=0$ otherwise.
	\end{example}
	
	\begin{example}[Square bipyramid: topological vs.\ graded Betti numbers]\label{example2}
		Let $\Delta$ be the simplicial complex on the vertex set $V=\{0,1,2,3,4,5\}$ obtained as follows.
		Declare the \emph{square base} on $\{1,2,3,4\}$ with edges
		\[
		\{1,2\},\{2,3\},\{3,4\},\{4,1\},
		\]
		and connect each base vertex to both apices $0$ and $5$:
		\[
		\{0,i\},\{5,i\}\quad\text{for }i\in\{1,2,3,4\}.
		\]
		The $2$-faces (facets) are the eight triangles
		\[
		\{0,1,2\},\{0,2,3\},\{0,3,4\},\{0,4,1\},\qquad
		\{5,1,2\},\{5,2,3\},\{5,3,4\},\{5,4,1\}.
		\]
		Combinatorially, this is the boundary of the square bipyramid (the octahedron). 
		
		Hence over any field $k$,
		\[
		\tilde H_0(\Delta;k)=0,\quad \tilde H_1(\Delta;k)=0,\quad \tilde H_2(\Delta;k)\cong k,
		\]
		so the (unreduced) topological Betti numbers are
		\[
		\beta_0=1,\qquad \beta_1=0,\qquad \beta_2=1.
		\]
		%		\paragraph{Stanley--Reisner ideal.}
		Let $S=k[x_0,\dots,x_5]$ be standard graded and write $I_\Delta\subset S$ for the Stanley--Reisner ideal of $\Delta$.
		The minimal nonfaces of $\Delta$ are precisely the missing edges
		\[
		\{0,5\},\ \{1,3\},\ \{2,4\},
		\]
		whence
		\[
		I_\Delta=(x_0x_5,\ x_1x_3,\ x_2x_4).
		\]
		% Minimal free resolution and graded Betti numbers.
		Because $I_\Delta$ is a complete intersection of three quadrics, the minimal free resolution of $S/I_\Delta$ over $S$ is the (pure) Koszul-type complex with shifts $2,2,2$:
		\[
		0 \longrightarrow 
		\underbrace{S(-6)}_{i=3}
		\longrightarrow 
		\underbrace{S(-4)^{3}}_{i=2}
		\longrightarrow 
		\underbrace{S(-2)^{3}}_{i=1}
		\longrightarrow 
		\underbrace{S}_{i=0}
		\longrightarrow S/I_\Delta \longrightarrow 0.
		\]
		Hence the graded Betti numbers are
		\[
		\beta_{0,0}=1,\quad \beta_{1,2}=3,\quad \beta_{2,4}=3,\quad \beta_{3,6}=1,
		\]
		and $\beta_{i,j}=0$ otherwise. Equivalently, the total Betti numbers are $(1,3,3,1)$.
	\end{example}
	
	
\subsection{Machine Learning Modeling}\label{sec:ML}

We train our regression models with 10-fold cross validation (CV) combined with Gradient Boosting Decision Trees (GBDT), implemented in Python via the \texttt{scikit-learn} library\,(v1.3.2). GBDT is valued for its resistance to overfitting, limited sensitivity to hyperparameter tuning, and straightforward deployment. It assembles many shallow decision trees generated from bootstrap samples of the training data and aggregates their outputs, so the ensemble corrects errors that any single tree might make. We supply the algorithm separately with PSRT-derived molecular descriptors and transformer-based descriptors. The hyperparameter settings are listed in \autoref{table:GBDT-parameters}.

\begin{table}[htb!]
	\small
	\centering
	\caption{Hyperparameters used in the gradient boosting regression model within a scikit-learn pipeline. StandardScaler is applied prior to model training.}
	\begin{tabular}{c c c c}
		\hline
		\textbf{No. of estimators} & \textbf{Max depth} & \textbf{Min. samples split} & \textbf{Learning rate} \\
		10,000 & 7 & 3 & 0.01 \\
		\hline
		\textbf{Max features} & \textbf{Subsample size} & \textbf{Random state} & \textbf{Standardization} \\
		Square root & 0.7 & Fixed (seeded) & Yes (StandardScaler) \\
		\hline
	\end{tabular}
	\label{table:GBDT-parameters}
\end{table}


\subsection{Evaluation Metrics}

We employ the Pearson correlation coefficient (\(R_p\)) for quantitative evaluation of the performance of our binding affinity prediction models, defined as: 
\begin{align*}
	\text{\(R_p\)}(\mathbf{x}, \mathbf{y}) = \frac{\sum_{m=1}^{M} (y_m^e - \bar{y}^e)(y_m^p - \bar{y}^p)}{\sqrt{\sum_{m=1}^{M} (y_m^e - \bar{y}^e)^2 \sum_{m=1}^{M} (y_m^p - \bar{y}^p)^2}},
\end{align*}
where \( y_m^e \) and \( y_m^p \) denote the experimental and predicted binding affinity values for the \( m \)-th sample, respectively, and \( \bar{y}^e \) and \( \bar{y}^p \) are their corresponding mean values.

We also report the root mean squared error (RMSE), which is computed as:
\begin{align*}
	\mathrm{RMSE} = \sqrt{\frac{1}{n} \sum_{m=1}^{M} (y_m^e - y_m^p)^2},
\end{align*}
where \( y_m^e \) and \( y_m^p \) represent the experimental and predicted binding affinity values for the \( m \)-th sample, respectively.

The above two metrics are employed to assess the performance of our machine learning models on all datasets. The original labels for these datasets are given as \( pK_d \) values, which are converted to binding free energies (in kcal/mol) by multiplying by a constant factor of 1.3633. Our models achieve reasonable RMSE values across all three datasets.


	\section{Conclusion}
	
	Commutative algebra provides a rigorous foundation for capturing combinatorial structure in biomolecular data. In this study, we introduce GBNL, a persistent commutative algebra framework for protein-DNA/RNA binding prediction.To our knowledge, this work presents the first graded-Betti application to large-data machine learning. DNA/RNA sequences are modeled as $k$-mer position sets, from which graded invariants of the associated Stanley-Reisner ideals are tracked over a positional filtration. In conjunction with transformer-based protein and nucleic acid embeddings, the algebraic features integrate nucleotide-level organization and protein-wide regularities into a concise, interpretable descriptor for predicting binding interactions.
	
	Beyond delivering precise affinity prediction, GBNL exhibits pronounced responsiveness to single-nucleotide mutations. For the N-China-F primer, a single nucleotide mutation yields a localized graded Betti shift that differentiates the mutated profile from the reference. A comprehensive mutation profiling of the same primer demonstrates that graded Betti signals capture subtle sequence perturbations beyond the resolution of standard Betti counts. Hence, GBNL upholds strong calibration and steady predictive behavior.. Furthermore, the GBNL achieves consistently strong performance by reaching a Pearson correlation of \(R_p = 0.69\) with RMSE \(= 1.81\)\,kcal/mol on S186, surpassing the SVSBI model~\cite{shen2023svsbi}, and maintaining comparable accuracy on S142 and S322. Collectively, these findings position GBNL as a broadly applicable and practical framework for biomolecular modeling.
	
	\section*{Data and Code availability}
	
	All codes needed to evaluate the conclusions in this study are available at \href{https://github.com/mzia-s/GBNL}{https://github.com/mzia-s/GBNL}. The three datasets can be found at \href{https://github.com/mzia-s/CAP}{https://github.com/mzia-s/CAP}. The SVSBI dataset can also be accessed at \href{https://github.com/WeilabMSU/SVS}{https://github.com/WeilabMSU/SVS}.
	
	\section*{Supporting Information}
	\href{./SI.pdf}{Supplementary Information} is available for supplementary tables. 
	
	\section*{Conflict of Interest}
	The authors declare no competing financial interests.
	
	
	\section*{Acknowledgments}
	This work was supported in part by NIH grant R35GM148196, NSF grant DMS-2052983, and MSU Research Foundation.    
	
	\clearpage
	
	
	\bibliographystyle{unsrt}
	\bibliography{refs}
	
	
	
\end{document}
