Introduction

Gene regulation is essential for controlling the timing and extent of gene expression in cells. It guides important processes from development to disease response.¹ Gene expression is tightly controlled by various DNA-binding proteins, including transcription factors (TFs) and epigenetic regulators.^2–4 TFs initiate gene expression by binding to specific DNA sequences, during which time a primary RNA transcript is synthesized from a gene’s DNA.⁵ On the other hand, epigenetic regulators control gene expression by binding to specific chromatin regions, which spread post-translational modifications that modulate 3D genome organization.⁶ Therefore, characterizing the DNA-interaction specificities of DNA binding proteins is critical for understanding the molecular mechanisms underlying many DNA-templated processes.⁷

To establish a comprehensive understanding of DNA binding processes, various experimental technologies have been developed and utilized.^8–16 These methods, such as Chromatin Immunoprecipitation followed by sequencing (ChIP-Seq),^8,9,17 Protein-binding microarray (PBM),^18–20 and Systematic Evolution of Ligands by Exponential Enrichment (SELEX),^15,16,21 have proven invaluable for measuring DNA-specific protein recognition. Nonetheless, due to the need for protein-specific antibodies,¹⁷ as well as the cost and intrinsic bias associated with these experiments,²² high-throughput measurement of large numbers of protein-DNA variants remains challenging.

Computational methods complement experimental efforts by providing the initial filter for assessing sequence-specific protein-DNA binding affinity. Numerous methods have emerged to enable predictions of binding sites and affinities of DNA-binding proteins. ^23–31 These methods often utilized machine-learning-based training to extract sequence preference information from DNA or protein by utilizing experimental high-throughput (HT) assays,^23–27,31 which rely on the availability and quality of experimental binding assays. Additionally, many approaches employ deep neural networks,^28–30 which could obscure the interpretation of interaction patterns governing protein-DNA binding specificities. Understanding these patterns, however, is crucial for elucidating the molecular mechanisms underlying various DNA-recognition processes, such as those seen in TFs.³²

Nowadays, over 5000 protein-DNA 3D structures, including TF-DNA complexes, have been published.^33,34 These data provide invaluable resources for understanding the physicochemical properties of protein-DNA binding patterns, extending beyond mere sequence information. Utilization of these data enables the training of a model for predicting the binding affinities or specificities of protein-DNA interactions. Emerging deep learning, ³⁵ probabilistic ³⁶ and chemistry-based models³⁷ are beginning to learn from protein–DNA cocrystal structures. Their predictions are biophysically meaningful and can be applied to understand molecular interactions at a mechanistic level, thus holding great potential to guide experimental design and synthetic biology applications. The very robustness of evolution^38–41 provides an opportunity to extract the sequence-structure relationships embedded in existing complexes. Guided by this principle, we can learn an interpretable binding energy landscape that governs the recognition processes of DNA-binding proteins.

Here, we introduce the Interpretable protein-DNA Energy Associative (IDEA) model, a predictive model that learns protein-DNA physicochemical interactions by fusing available biophysical structures and their associated sequences into an optimized energy model (Figure 1). We show that the model can be used to accurately predict the sequence-specific DNA binding affinities of DNA-binding proteins and is transferrable across the same protein superfamily. Moreover, the model can be enhanced by incorporating experimental binding data and can be generalized to enable base-pair resolution predictions of genomic DNA-binding sites. Notably, IDEA learns a family-specific interaction matrix that quantifies energetic interactions between each amino acid and nucleotide, allowing for a direct interpretation of the “molecular grammar” governing sequence-specific protein-DNA binding affinities. This interpretable energy model is further integrated into a simulation framework, enabling mechanistic studies of various biomolecular functions involving protein-DNA dynamics.

Figure 1:

Results

Predictive Protein-DNA Energy Model at Residue Resolution

IDEA is a coarse-grained biophysical model at the residue resolution for investigating protein-DNA binding interactions (Figure 1). It integrates both structures and corresponding sequences of known protein-DNA complexes to learn an interpretable energy model based on the interacting amino acids and nucleotides at the protein-DNA binding interface. The model is trained using available protein-DNA complexes curated from existing databases. ^33,42 Unlike existing deep-learning-based protein-DNA binding prediction models, IDEA aims to learn a physicochemical-based energy model that quantitatively characterizes sequence-specific interactions between amino acids and nucleotides, thereby interpreting the “molecular grammar” driving the binding energetics of protein-DNA interactions. The optimized energy model can be used to predict the binding affinity of any given protein-DNA pair based on its structures and sequences. Additionally, it enables the prediction of genomic DNA binding sites by a given protein, such as a transcription factor. Finally, the learned energy model can be incorporated into a simulation framework to study the dynamics of DNA-binding processes, revealing mechanistic insights into various DNA-templated processes. Further details of the optimization protocol are provided in the Methods Section Energy Model Optimization.

IDEA Accurately Predicts Sequence-Specific Protein-DNA Binding Affinity

We first examine the predictive accuracy of IDEA by comparing its predicted TF-DNA binding affinities with experimental measurements.^13,14 We focused on the MAX TF, a basic Helix-loop-helix (bHLH) TF with the most comprehensive available experimental binding data. The binding affinity of MAX to various DNA sequences has been quantified by multiple experimental platforms, including different SELEX variants^15,16,21 and microfluidic-based MITOMI assay.^13,43 Among them, MITOMI quantitatively measured the binding affinities of MAX TF to a comprehensive set of 255 DNA sequences with mutations in the enhancer box (E-box) motif, the consensus MAX DNA binding sequence. This dataset serves as a reference for us to benchmark our model predictions. Our de novo prediction, based on one MAX crystal structure (PDB ID: 1HLO), correlates well with the experimental values (Pearson correlation coefficient 0.67) and correctly ranks the binding affinities of those DNA sequences (Spearman’s rank correlation coefficient 0.65, Figure 2A). Including additional human-associated MAX-DNA complex structures and their associated sequences in the model training slightly improves the predictions (Pearson correlation coefficient 0.68, Spearman’s rank correlation coefficient 0.65, Figure S1), albeit not significantly, likely due to the high structural similarity of the protein-DNA interfaces across all MAX proteins. Prior works have proven it instrumental in incorporating experimental binding data to improve predicted protein-nucleic acid binding affinities.^24,27,44 Inspired by these approaches, we developed an enhanced predictive protocol that integrates additional experimental binding data if available (see the Methods Section Enhanced Modeling Prediction with SELEX Data). Encouragingly, when training the model with the SELEX-seq data²⁷ of MAX TF, IDEA showed a significant improvement in reproducing the MITOMI measurement (Pearson correlation coefficient 0.79, Spearman’s rank correlation coefficient 0.76, Figure S2).

Figure 2:

IDEA Decodes the Molecular Grammar Governing Protein-DNA Interactions

IDEA protocol learns the physicochemical interactions that determine protein-DNA interactions by utilizing the sequence-structure relationship embedded in the protein-DNA experimental structures. Such a physicochemical interaction pattern can be interpreted directly from the learned energy model. To illustrate this, we focused on the energy model learned from the MAX-DNA complexes (Figure 2B). Notably, DNA deoxycytidine (DC) exhibited strong interactions with protein arginine (R), consistent with the fact that the E-box region (CACGTG) frequently attracts the positively charged residues of bHLH TFs.⁴⁵ Importantly, arginine was in close contact with deoxycytidine in the crystal structure (Figure 2C), thus consistent with the strong DC-R interactions shown in the learned energy model. Upon integrating the SELEX data into our model training, we found that the improved energy model shows additional unfavorable interactions between protein glutamic acid (E) and deoxycytidine, consistent with their negative charges (Figure S3). Thus, including more experimental data can boost IDEA’s predictive accuracy by refining the amino-acid-nucleotide interacting energy model to better align with physical principles.

IDEA Generalizes across Various Protein Families

To examine IDEA’s predictive accuracy across different DNA-binding protein families, we applied it to calculate protein-DNA binding affinities using a comprehensive HT-SELEX dataset.²⁶ We focused on evaluating the capability of IDEA to distinguish strong binders from weak binders for each protein with an experimentally determined structure. We calculated the probability density distribution of the top and bottom binders identified in the SELEX experiment. A well-separated distribution indicates the successful identification of strong binders by IDEA (Figure 2D and S4). Receiver Operating Characteristic (ROC) and Precision-Recall (PR) analyses were performed to calculate the Area Under the Curve (AUC) and the precision-recall AUC (PRAUC) scores for these predictions. Further details are provided in the Methods Section Evaluation of IDEA Prediction Using HT-SELEX Data. Our analysis shows that IDEA successfully differentiates strong from weak binders for 80% of the 22 proteins across 12 protein families, achieving AUC and balanced PRAUC scores greater than 0.5 (Figure 2E and S5). To benchmark IDEA’s performance against other leading methods, we compared its predictions with several popular models, including the sequence-based predictive models ProBound³¹ and DeepBind,²⁴ the family-based energy model rCLAMPS,³⁶ and the knowledge-based energy model DBD-Hunter.⁴⁶ IDEA demonstrates performance comparable to these state-of-the-art approaches, and incorporating sequence features further improves its prediction accuracy (Figure S6, Table S1, and Table S2). We also performed 10-fold cross-validation on the binding affinities of protein–DNA pairs in this dataset and found that IDEA outperforms a recent regression model that considers the shape of DNA with different sequences²⁶ (Figure S7). Details are provided in the Appendix Section: Comparison of IDEA predictive performance Using HT-SELEX data.

We also applied IDEA to predict the binding affinity of additional TFs with available MITOMI measurements.⁴³ For PHO4, another bHLH transcription factor, IDEA’s predictions, trained on the only available protein-DNA structure (PDB ID: 1A0A), correlate well with experimental measurements, showing a Pearson correlation coefficient of 0.60 and a Spearman’s rank correlation coefficient of 0.63 (Figures S8A, B). We further evaluated IDEA’s predictive performance for the zinc-finger protein Zif268 (see the Methods Section Structural Modeling of Protein and DNA). Due to the limited experimental data for all possible DNA sequence combinations, we focused on testing IDEA’s predictions on point-mutated DNA sequences. Predictions on point-mutated sequences are known to be a challenging task due to their minor deviations from the wild-type sequence. Despite this, IDEA achieves accurate predictions, with a Pearson correlation coefficient of 0.57 and a Spearman’s rank correlation coefficient of 0.60 (Figures S8C, D). We further expanded the training dataset to include all available Zif268 structures and their associated sequences from the same CATH zinc finger superfamily (CATH ID: 3.30.160.60). Incorporating these additional training data further improves the predictive accuracy (Pearson correlation coefficient 0.63; Spearman’s rank correlation coefficient 0.60, Figure S9).

IDEA Demonstrates Transferability across Proteins in the Same CATH Superfamily

Since IDEA relies on the sequence-structure relationship of given protein-DNA complexes to reach predictive accuracy, we inquired whether the trained energy model from one protein-DNA complex could be generalized to predict the sequence-specific binding affinities of other complexes. To test this, we assessed the transferability of IDEA predictions across all 11 structurally available protein-DNA complexes within the MAX TF-associated CATH super-family (CATH ID: 4.10.280.10, Helix-loop-helix DNA-binding domain). We trained IDEA based on each of these 11 complexes and then used the trained model to predict the MAX-based MITOMI binding affinity. Our results show that IDEA generally makes correct predictions of the binding affinity when trained on proteins that are homologous to MAX, with Pearson and Spearman Correlation coefficients larger than 0.5 (Figure 3A and Figure S10).

Figure 3:

The transferability of IDEA within the same CATH superfamily can be understood from the similarities in protein-DNA binding interfaces, which determine similar learned energy models. For example, the PHO4 protein (PDB ID: 1A0A) shares a highly similar DNA-binding interface with the MAX protein (PDB ID: 1HLO) (Figure 3B), despite sharing only a 33.41% probability of being homologous. Consequently, the energy model derived from the PHO4-DNA complex (Figure 3C) exhibits a similar amino-acid-nucleotide interactive pattern as that learned from the MAX-DNA complex (Figure 2B). To further evaluate the similarity between the learned energy models and their connection to protein families, we performed principal component analysis (PCA) on the normalized energy models across 24 proteins from 12 protein families.²⁶ Our analysis (Figure S11) reveals that most of the energy models from the same protein family fall within the same cluster, while those from different protein families exhibit distinct patterns. Moreover, the relative distance between energy models in PCA space reflects the degree of transferability between them. For example, PHO4 (PDB ID: 1A0A) is positioned close to MAX (PDB ID: 1HLO), whereas USF1 (PDB ID: 1AN4) and TCF4 (PDB ID: 6OD3) are farther away. This is consistent with the results in Figure 3A, where the energy model trained on PHO4 has better transferability than those trained on USF1 or TCF4.

Identification of Genomic Protein-DNA Binding Sites

The genomic locations of DNA binding sites are causally related to major cellular processes. ¹⁷ Although multiple techniques have been developed to enable high-throughput mapping of genomic protein-DNA binding locations, such as ChIP–seq, ^9,17,49 DAP-seq,¹⁰ and FAIRE-seq,⁵⁰ challenges remain for precisely pinpointing protein-binding sites at a base-pair resolution. Furthermore, the applicability of these techniques for different DNA-binding proteins is restricted by the quality of antibody designs.¹⁷ Therefore, a predictive computational framework would significantly reduce the costs and accelerate the identification of genomic binding sites of DNA-binding proteins.

We incorporated IDEA into a protocol to predict the genomic protein-DNA binding sites at the base-pair resolution, given a DNA-binding protein and a genomic DNA sequence. To evaluate the predictive accuracy of this protocol, we utilized publicly available ChIP-seq data for MAX TF binding in GM12878 lymphoblastoid cell lines.⁵¹ As the experimental measurements were conducted at a 420 base pairs resolution, we averaged our modeling prediction over a window spanning 500 base pairs. Since the experimental signals are sparsely distributed across the genome, we focused our prediction on the 1 Mb region of Chromosome 1, which has the densest and most reliable ChIP-seq signals (156000 kb – 157000 kb, representative window shown in Figure 4A top). In our model, lower predicted binding free energy corresponds to stronger protein-DNA binding affinity and thus a higher probability of a binding site. These predicted binding free energies were further normalized against those calculated from the randomized decoy sequences (see the Methods Section Processing of ChIP-seq Data for details). The identified MAX-binding sites are marked in red, with a normalized Z score < −0.75. Our prediction successfully identifies the majority of the MAX binding sites, achieving an AUC score of 0.81 (Figure 4B). We also tested IDEA’s predictive performance on two additional typical DNA-binding proteins and successfully identified their genomic binding sites (Figure S12). Therefore, IDEA provides an accurate initial prediction of genomic recognition sites of DNA-binding proteins based only on DNA sequences, facilitating more focused experimental efforts. Furthermore, IDEA holds the potential to identify binding sites beyond the experimental resolution and can aid in the interpretation of other sequencing techniques to uncover additional factors that influence genome-wide DNA binding, such as chromatin accessibility.^52–54

Figure 4:

Integration of IDEA into a Residue-Resolution Simulation Model Captures Protein-DNA Binding Dynamics

The trained residue-level model can be incorporated into a simulation model, thus enabling investigations into dynamic interactions between protein and DNA for mechanistic studies. Coarse-grained protein-nucleic acid models have shown strong advantages in investigating large biomolecular systems.^55–59 Notably, the structure-based protein model⁶⁰ and 3-Site-Per-Nucleotide (3SPN) DNA model⁶¹ have been developed at the residue resolution to enable efficient sampling of biomolecular systems involving protein and DNA. A combination of these two models has proved instrumental in quantitatively studying multiple large chromatin systems,⁶² including nucleosomal DNA unwrapping,^63–65 interactions between nucleosomes,^66,67 large chromatin organizations,^67–69 and its modulation by chromatin factors.^70–72 Although both the protein and DNA force fields have been systematically optimized,^61,73–76 the non-bonded interactions between protein and DNA were primarily guided by Debye-Hückel treatment of electrostatic interactions, which did not consider sequencespecific interactions between individual protein amino acids and DNA nucleotides.

To refine this simulation model by including protein-DNA interaction specificity, we incorporated the IDEA-optimized protein-DNA energy model into the coarse-grained simulation model as a short-range Van der Waals (vdW) interaction in the form of a Tanh function (See Methods Section Sequence Specific Protein-DNA Simulation Model for modeling details). This refined simulation model is used to simulate the dynamic interactions between DNA-binding proteins and their target DNA sequences. As a demonstration of the working principle of this modeling approach, we selected 9 TF-DNA complexes whose binding affinities were measured by the same experimental group.⁷⁷ We applied IDEA to train an energy model based on these protein-DNA complexes, which yielded a modest correlation between the modeling-predicted binding affinity and experimental measurements (Pearson correlation coefficient 0.46, Spearman’s rank correlation coefficient 0.55, Figure 5A). However, upon incorporating the learned energy model into our simulation model and computing protein-DNA binding free energies through umbrella sampling, ⁷⁸ we observed a significant improvement in predictive accuracy. The predicted binding free energies⁷⁹ (Figure S13) from our simulations correlate strongly with experimental measurements (Pearson correlation coefficient 0.81; Spearman’s rank correlation coefficient 0.90) and near-perfectly reproduce the experimental binding affinities (Figure 5B). This prediction also represents a substantial improvement over the previous model, which depicts protein-DNA interactions through non-sequence-specific homogeneous electrostatic interactions (Pearson correlation coefficient 0.60, Spearman’s rank correlation coefficient 0.70; Figure S14). Notably, due to an undetermined prefactor in the IDEA training, it can only provide relative binding free energies in reduced units for different protein-DNA complexes. In contrast, the IDEA-incorporated simulation model can predict the absolute binding free energies with physical units.

Figure 5:

Discussion

The protein-DNA interaction landscape has evolved to facilitate precise targeting of proteins towards their functional binding sites, which underlie essential processes in controlling gene expression. These interaction specifics are determined by physicochemical interactions between amino acids and nucleotides. By integrating sequences and structural data from available protein-DNA complexes into an interaction matrix, we introduce IDEA, a data-driven method that optimizes a system-specific energy model. This model enables high-throughput in silico predictions of protein-DNA binding affinities and can be scaled up to predict genomic binding sites of DNA-binding proteins, such as TFs. IDEA achieves accurate de novo predictions using only protein-DNA complex structures and their associated sequences, but its accuracy can be further enhanced by incorporating available experimental data from other binding assay measurements, such as the SELEX data,^15,16,21 achieving accuracy comparable or better than state-of-the-art methods (Figures S2 and S7, Table S1 and S2). Despite significant progress in genome-wide sequencing techniques, ^9,10,17,50 determining sequence-specific binding affinities of DNA-binding biomolecules remains time-consuming and expensive. Therefore, IDEA presents a cost-effective alternative for generating the initial predictions before pursuing further experimental refinement.

A key advantage of IDEA is its incorporation of both structural and sequence information, which greatly reduces the demand for extensive training data. Prior efforts have applied deep learning techniques to predict the interactions between proteins and DNA based on their sequences.^24,25,27 These approaches usually require a large amount of sequencing training data. By leveraging the sequence-structure relationship inherent in protein-DNA complex structures, IDEA achieves accurate de novo prediction using only a few protein-DNA complexes within the same protein superfamily.

In addition, despite the complicated in vivo environment, IDEA enables predictions of genomic binding sites and shows good agreement with the ChIP-seq measurements (Figure 4, S12). Moreover, IDEA features rapid prediction, typically requiring 1-3 days to predict base-pair resolution genomic binding sites of a 1 Mb DNA region using a single CPU. Furthermore, IDEA’s prediction can be trivially parallelized, making it possible to perform genome-wide predictions of DNA recognition sites by a given protein within weeks.

Another highlight of IDEA is its ability to present an interpretable, family-specific amino acid-nucleotide interaction energy model for given protein-DNA complexes. The optimized IDEA energy model can not only predict sequence-specific binding affinities of protein-DNA pairs but also provide a residue-specific interaction matrix that dictates the preferences of amino acid-nucleotide interactions within specific protein families (Figure S11). This interpretable energy matrix would facilitate the discovery of sequence binding motifs for target DNA-binding proteins, complementing both sequence-based^80–82 and structure-based approaches.^35–37,83 Additionally, we integrated this physicochemical-based energy model into a simulation framework, thereby improving the characterization of protein-DNA binding dynamics. IDEA-based simulation enables the investigation into dynamic interactions between various proteins and DNA, facilitating molecular-level understanding of the physical mechanisms underlying many DNA-binding processes, such as transcription, epigenetic regulations, and their modulation by sequence variations, such as single-nucleotide polymorphisms (SNPs).^84,85

Since the IDEA-optimized energy model serves as a “molecular grammar” for guiding protein-DNA interaction specifics, it can be expanded to include additional variants of amino acids and nucleotides. With recent advancements in sequencing and structural characterization,^86,87 as well as deep-learning-based structural predictions,^88–90 thousands of structures and sequences have been solved for these variants. The IDEA procedure can be repeated for these variant-included structures to expand the optimized energy model by including modifications such as methylated DNA nucleotides ⁹¹ and post-translationally modified amino acids.⁹² Such strategies will facilitate studies on the functional relevance of epigenetic variants, such as those caused by exposure to environmental hazards,⁹³ and their structural impact on the human genome.^94–96

Although IDEA has proven successful in many examples, it can be improved in several aspects. The model currently assumes the training and testing sequences share the same protein-DNA structure. While double-stranded DNA is generally rigid, recent studies have shown that sequence-dependent DNA shape contributes to their binding specificity.^25,35,97 To improve predictive accuracy, one could incorporate predicted DNA shapes or structures into the IDEA training protocol. In addition, the model is residue-based and evaluates the binding free energy as the additive sum of contributions from individual amino-acidnucleotide contacts. This assumption does not account for cooperative effects that may arise from multiple nucleotide changes. A potential refinement could utilize a finer-grained model that includes more atom types within contacting residues and employs a many-body potential to account for such cooperative effects.

Methods

Structural Modeling of Protein and DNA

We utilized a coarse-grained framework to extract structural information of protein and DNA at the residue level (Figure 1). Specifically, each protein amino acid was represented by the coordinates of the Cα atom, and each DNA nucleotide was represented by either the P atom in the phosphate group or the C5 atom on the nucleobase. The choice between these two DNA atom types (C5 or P) depended on their distances to surrounding Cα atoms. We hypothesized that DNA atoms closer to Cα atoms would provide more meaningful information on protein-DNA interactions. Therefore, we chose the DNA atom type (C5 or P) with the largest sum of occurrences within 8°A, 9°A, and 10°A distances from their surrounding Cα atoms (Table S3). This selection rule was applied consistently across all complexes in the training dataset to maintain uniform criteria.

Energy Model Optimization

Training Protocol

IDEA is a residue-based, data-driven biophysical model that predicts sequence-specific protein-DNA binding affinities. The model outputs relative binding free energies between protein and DNA, which can be converted to the relative binding affinity (represented by the dissociation constant K_d) using:⁷

The model integrates both the sequence and structural information from experimentally determined protein-DNA structural interface into an interaction matrix for model training. These sequences, located at the protein-DNA structural interface (defined as those residues whose representative atoms have a distance ≤ 8°A, highlighted in blue in Figure 1 Left), are hypothesized to have evolved to favor strong amino acid-nucleotide interactions. We collected the sequences from both the strong binders and decoy binders (considered weak binders) for model parameterization. Synthetic decoy binders were generated by randomizing either the DNA (1000 sequences) or protein sequences (10,000 sequences) of the strong binders (see Figure S15 for an assessment of the model’s robustness in transferability and generalizability across different numbers of decoys).

IDEA optimizes an amino-acid-nucleotide energy matrix, γ(a, n), by maximizing the gap in binding energies between the strong and decoy binders. Similar strategies have been successfully applied in protein folding and protein-protein binding.^{38,40,56,98,99} Specifically, the solvent-averaged binding free energies for strong binders (E_strong) and their corresponding decoy weak binders (E_decoy) are calculated using the following expression:

Here, γ(a, n) is a 20-by-4 residue-type-specific energy matrix encoding residue-type-specific interaction energies between amino acid (a) and nucleotide (n). Θ_i,j(a_i, n_j) utilizes a switching function that captures an effective interaction range between each amino acid-nucleotide pair within the protein-DNA complex:

here, r_min = −8°A and r_max = 8°A. This defines two residues to be “in contact” if their distance is less than 8°A. κ (here taken 0.7) is a scaling factor that modulates the steepness of the hyperbolic tangent function.

In practice, the γ matrix was optimized to maximize δE/ΔE, where δE = ⟨E_decoy⟩ − ⟨E_strong⟩ represents the average energy gap between strong and decoy binders, and ΔE is the standard deviation of decoy binding energies. Mathematically,

and

with

Here, summarizes the total number of contacts between each type of amino acid and nucleotide in the given protein-DNA complexes. The vector A(a, n) quantifies the difference in interaction frequency between the strong and decoy binders, with a dimension of 1 × 300. The covariance matrix B(a, n) is a 300 × 300 matrix that quantifies the propensity of co-occurring residue type interactions in the decoy binders. Maximizing the function δE/ΔE with respect to γ(a, n) corresponds to maximizing , which can be achieved by maximizing the functional objective R(γ):

where λ is a Lagrange multiplier. Setting the derivative equal to zero leads to γ ∝ B⁻¹A. A filtering procedure is further applied to reduce the noise arising from the finite sampling of certain types of interactions: The matrix B is first diagonalized such that B⁻¹ = P Λ⁻¹P ⁻¹, where P is the matrix of eigenvectors of B, and Λ is a diagonal matrix composed of B’s eigenvalues. We retained the principal N modes of B, which are ordered by the descending magnitude of their eigenvalues and replaced the remaining 300 − N eigenvalues in Λ with the N th eigenvalue of B. In all optimizations reported in this work, N = 70 was generally a good choice to maximize the utilization of non-zero eigenvalues in the matrix B. For visualization purposes, the vector γ is reshaped into a 20-by-4 matrix, as shown in Figure 2B.

To predict the binding free energy for a target protein-DNA pair, we substitute the target sequence into a known protein-DNA complex structure and recalculate the target ϕ matrix, which will be combined with the trained γ matrix to compute the target protein-DNA binding free energies using Eq. 2. Due to the presence of an undetermined scaling factor after the model optimization, the predicted free energies are presented in reduced units. Despite this, the relative free energies can be used to accurately rank binding affinities of target protein-DNA pairs. We utilize this computed solvent-averaged free energy to approximate the protein-DNA binding free energy in Equation 1. When the optimized energy is integrated into a simulation model, the IDEA-based simulations can predict the absolute binding free energies in physical units (Figure 5C).

Treatment of Complementary DNA Sequences

To replace the DNA sequence in the protein-DNA complex structure with a target sequence, IDEA uses sequence identity to determine whether the target sequence belongs to the forward or reverse strand of the DNA in the protein-DNA structure. The more similar strand is selected and replaced with the target sequence. As the orientations of test sequences are specified from 5’ to 3’ in all datasets used in this study (e.g., processed MITOMI, HT-SELEX, and ChIP-seq data), this approach ensures that the target sequences are replaced and evaluated correctly. In cases where sequence orientation is not provided (though this was not an issue in this study), we recommend replacing both the forward and reverse strands with the target sequence separately and evaluating the corresponding protein–DNA binding free energies. Since strong binders are likely to dominate the experimental signals, the higher predicted binding affinity, with stronger binding free energies, should be taken as the model’s final prediction.

Enhanced Modeling Prediction with SELEX Data

When additional experimental binding data, such as SELEX,²⁷ is available, we extended IDEA to incorporate this data for enhancing the model’s predictive capabilities. SELEX data provides binding affinities between given transcription factors and DNA sequences, which can be converted into dimensionless binding free energy using Eq. 1, with K_d being the normalized affinity generated from the SELEX package. We then calculated the ϕ value by utilizing the protein-DNA complex structure, replacing the native DNA sequence with the SELEX-provided sequences. These ϕ values can be used together with the converted ΔΔG values to compute the γ energy model based on the equation 2:

here γ is the enhanced protein-DNA energy model represented as a 300 × 1 vector. We estimated γ by applying ridge regression to the linear system in Equation 9, using the previously computed ϕ matrix and the SELEX-derived ΔΔG values. A regularization parameter of α = 0.01 was consistently applied throughout all computations to prevent overfitting.

For predicting the MAX-DNA binding specificity, we utilized the Human MAX SELEX data deposited in the European Nucleotide Archive database (accession number PRJEB25690).^27,100 The R/Bioconductor package SELEX version 1.30.0¹⁰¹ was used to determine the observed R1 count for all 10-mers. Among them, ACCACGTGGT is the motif with the highest frequency, which aligns with the DNA sequence in the MAX-DNA PDB structure (PDB ID: 1HLO), whose full DNA sequence is CACCACGTGGT. Consequently, we used this motif as the reference to align the remaining SELEX-measured 10-mer sequences. A total of 255 10-mer sequences, corresponding to the DNA variants measured from the MIT-OMI experiment,¹³ were selected. We threaded the PDB structure with these sequences to construct the ϕ and utilized the SELEX-converted binding free energies to compute the γ energy model.

Evaluation of IDEA Predictions Using HT-SELEX data

The processed HT-SELEX data used in this study are from,²⁶ which contains processed DNA binding sequence motifs and their normalized copy number detected from the HT-SELEX experiment, referred to as M-word scores. A higher M-word score corresponds to a stronger binding sequence motif detected in the HT-SELEX experiments. Among all the proteins in the data, 23 have experimentally determined protein-DNA complex structures. We predicted the binding free energies of all the processed sequences of these proteins using our protocol. For evaluating IDEA predictions, we classified motifs with M-word scores above 0.9 as strong binders (label 1) and those with M-word scores below 0.3 as weak binders (label 0). In cases where no sequences had an M-word score below 0.3, we used alternative cutoffs of 0.4 or 0.5 to ensure that at least three weak-binding sequences were included. These labeled strong and weak binders were then used as the ground truth for evaluating IDEA’s predictions. We assessed IDEA’s predictive accuracy by calculating the probability density distributions of the predicted binding free energies (Figure S4). A well-separated distribution between strong and weak binders indicates a successful prediction. Additionally, we calculated the AUC score and balanced PRAUC score based on the Receiver Operating Characteristic (ROC) and Precision-Recall (PR) analyses of these probability distributions. More details about the evaluation metrics and comparisons with popular existing models are provided in the Appendix Section: Comparison of IDEA predictive performance Using HT-SELEX data

Selection of CATH Superfamily for Testing Model Transferability

The MAX protein is characterized by a Helix-loop-helix DNA-binding domain and belongs to the CATH superfamily 4.10.280.10,¹⁰² which includes a total of 11 protein-DNA structures with the E-box sequence (CACGTG), similar to the MAX-binding DNA. We hypothesized that proteins within the same superfamily exhibit common residue-interaction patterns due to their structural and evolutionary similarities. To test the transferability of our model, we conducted individual training on each one of those 11 protein-DNA complex structures and their associated sequences, predicting the MAX-DNA binding affinity in the MITOMI dataset.¹³ To examine the connection between the predictive outcome and the probability that the training protein is homologous to MAX, HHpred¹⁰³ was utilized to search for MAX’s homologous proteins, using the sequence of 1HLO as the query (Figure 3A).

Processing of ChIP-seq Data

To evaluate the predictive capabilities of our protocol on specific chromosomal segments, we selected ChIP-seq data for three transcription factors—MAX, EGR1, and CTCF. The data were sourced from the ENCODE project, with accession numbers ENCFF361EVH (MAX), ENCSR026GSW (EGR1), and ENCFF960ZGP (CTCF).

For testing, we chose genomic regions with the densest ChIP-seq signals. For MAX and CTCF, we focused on the human GM12878 chromosome 1 (GRCh38) region between 156,000,000 bp and 157,000,000 bp, as this 1 Mb segment contains the densest signals for both transcription factors, with 33 MAX binding sites and 70 CTCF binding sites. For EGR1, we used the HepG2 cell line, focusing on the GRCh38 chr8 region between 143,000,000 bp and 144,000,000 bp, which contains 97 peaks.

To predict binding sites, we scanned the corresponding DNA regions base by base using the DNA sequence closest to each protein in the relevant protein-DNA complex structures. For MAX, we used the 11-bp sequence from the PDB structure 1HLO, generating 1, 000, 000− 11 + 1 = 999, 990 sequences. For EGR1, we used the 10-bp sequence from PDB 1AAY, producing 1, 000, 000 − 10 + 1 = 999, 991 sequences. For CTCF, no single PDB structure covers the entire protein,¹⁰⁴ so two PDB structures were used: 8SSS (ZF1-ZF7, 23-bp) and 8SSQ (ZnFs 3-11, 35-bp), generating 1, 000, 000 −23 + 1 = 999, 978 and 1, 000, 000 −35 + 1 = 999, 966 sequences, respectively.

Predicted binding free energies for each transcription factor were normalized into Z-scores using randomized decoy sequences:

where ⟨ ⟩ is the average over all decoy sequences, and SD is the standard deviation. Given the evolutionary preference for strong binding, most Z-scores were negative, indicating stronger binding. ROC curves were then constructed to assess the predictive performance, with lower predicted Z-scores representing stronger binding affinity (Figure 4, S11).

Sequence Specific Protein-DNA Simulation Model

We utilized a previously developed protein and DNA residue-level coarse-grained model for simulating protein-DNA binding interactions. The protein was modeled using the Cα-based structure-based model,⁶⁰ with each amino acid represented by a simulation bead located at the Cα site. The DNA was modeled with the 3-Site-Per-Nucleotide (3SPN) DNA model, ⁶¹ with each DNA nucleotide modeled by three coarse-grained sites located at the centers of mass of the phosphate, sugar, and base sites. Both these two models and their combination have been validated to reproduce multiple experimental measurements.^{60,61,64,66–69,71,72,75,105}

Detailed descriptions of this protein-DNA force field have been documented in previous studies.^63,68,69,72 Specifically, we followed Ref. 61 for simulating DNA-DNA interactions, and Ref. 60 for simulating the protein-protein interactions. The protein-DNA interactions were modeled with a long-range electrostatic interaction and a short-range sequence-specific interaction. The electrostatic interactions were simulated with the Debye-Hückel treatment to capture the screening effect at different ionic concentrations:

where ɛ = 78.0 is the dielectric constant of the bulk solvent. q_i and q_j correspond to the charges of the two particles. ɛ₀ = 1 is the vacuum electric permittivity, and l_D is the Debye-Hückel screening length.

The excluded volume effect was modeled using a hard wall potential with an apparent radius of σ = 4.0°A:

An additional sequence-specific protein-DNA interaction was modeled with a Tanh function between the phosphate site of the DNA and the amino acid, with the form:

where r₀ = 8°A and η = 0.7°A⁻¹. W is a 4 × 20 matrix that characterizes the interaction energy between each type of amino acid and DNA nucleotide. Here, We used the normalized IDEA-learned energy model based on 9 protein-DNA complexes (Figure 5A). The IDEA-learned energy model was normalized by the sum of the absolute value of each term of the energy model. Detailed parameters are provided in Table S4. Additional details on the protein-DNA simulation are provided in the Supplementary Materials.

Supplementary Figures

Figure S1:

Figure S2:

Figure S3:

Figure S4:

Figure S5:

Figure S6:

Figure S7:

Figure S8:

Figure S9:

Figure S10:

Figure S11:

Figure S12:

Figure S13:

Figure S14:

Figure S15:

Supplementary Tables

Table S1:

Table S2:

Table S3:

Table S4:

Data availability

The implementation of the IDEA model, along with training and prediction examples, is available on our GitHub repository and in the Supplementary data.

Acknowledgements

This work was supported by the startup funding from North Carolina State University. Additionally, this research received partial funding from the NC State Genetics and Genomics Academy and Comparative Medicine Institute. We appreciate Dr. Keith Weninger, Dr. Faruck Morcos, and Dr. Qin Zhou for their insightful discussions and Dr. Sebastian Maerkl for guiding us to the MITOMI experimental data.

Additional files

Supplementary Information

Appendix

Details of Molecular Dynamics Simulation

The simulations were carried out using the LAMMPS software¹⁰⁶ with the implementation of the protein-DNA force field (Table S4). We performed umbrella-sampling simulations ⁷⁸ for comprehensive sampling of the binding and unbinding events to compute the protein-DNA binding affinities. To evaluate the model’s performance, we selected 9 protein-DNA complexes with binding affinity measured by the same experimental lab.^77,107–110 Initial simulation configurations were prepared based on the crystal structures of the protein-DNA complexes. The PDB IDs of these complexes are provided in Figs 5C, S13, and S14.

The simulations were conducted under the same experimental conditions with 100 mM monovalent ions. We used a spring constant of 0.02 kcal/mol/°A² based on the center of mass (COM) distances between the protein and DNA. The centers of the umbrella windows were placed on a grid ranging from 0.0 to 140.0°A, with a grid size of 10.0°A, and each umbrella simulation lasted for 22 million steps, with a time step of 2.0 fs. We simulated 6 copies of the simulation with different random seeds to ensure sampling convergence. When reweighting the simulation to compute the free energy⁷⁹ (Figs. S13 and S14), we excluded the first 2 million steps of each simulation trajectory. To calculate the error bars from our simulation, we divided each simulation trajectory into three equal-length blocks and independently calculated the reweighted free energies from them. The error bar was then estimated as the standard deviation of the three free energy profiles.

Comparison of IDEA predictive performance Using HT-SELEX data

To benchmark the performance of IDEA against state-of-the-art protein-DNA predictive models, we evaluated its ability to recognize strong binders with the HT-SELEX datasets across 22 proteins from 12 families.²⁶ Specifically, we compare IDEA with two widely used sequence-based models: ProBound³¹ and DeepBind.²⁴ ProBound has demonstrated superior performance over many other predictive protein-DNA models, including JASPAR 2018, ¹¹¹ HOCOMOCO,¹¹² Jolma et al.,¹¹³ and DeepSELEX.¹¹⁴ To use ProBound, we retrieved the trained binding model for each protein from motifcentral.org and used the GitHub implementation of ProBoundTools to infer the binding scores between protein and target DNA sequences. Except for POU3F1, binding models are available for all proteins. Therefore, we excluded POU3F1 and evaluated the protein-DNA binding affinities for the remaining 21 proteins. To use DeepBind, sequence-specific binding affinities were predicted directly with its web server. The Area Under the Curve (AUC) and the Precision-Recall AUC (PRAUC) scores were used as metrics for comparison. An AUC score of 1.0 indicates a perfect separation between the strong– and weak-binder distributions, while an AUC score of 0.5 indicates no separation. Because there is a significant imbalance in the number of strong and weak binders from the experimental data,²⁶ where the strong binders are far fewer than the weak binders, we reweighted the samples to achieve a balanced evaluation, using 0.5 as the base-line for randomized prediction.¹¹⁵ As summarized in Figure S6, Table S1, and Table S2, IDEA ranked second among the three predictive models. In order to assess the performance of IDEA when augmented with additional protein-DNA binding data, we augmented IDEA using randomly selected half of the HT-SELEX data (see the Methods Section Enhanced Modeling Prediction with SELEX Data). The augmented IDEA model achieved the best performance among all the models.

In addition, we compared the performance of IDEA with both general and family-specific knowledge-based energy models. First, we incorporated a knowledge-based generic protein-DNA energy model (DBD-Hunter) learned from the protein-DNA database, reported by Skoinick and coworkers,⁴⁶ into our prediction protocol. This model assigns interaction energies to different functional groups within each DNA nucleotide, including phosphate (PP), sugar (SU), pyrimidine (PY), and imidazole (IM) groups. For our comparison, we averaged the energy contributions of these groups within each nucleotide and replaced the IDEA-learned energy model with the DBD-Hunter model to assess its ability to differentiate strong binders from weak binders in the HT-SELEX dataset. ²⁶ Additionally, we compared IDEA with rCLAMPS, a family-specific energy model developed to predict protein-DNA binding specificity in the C2H2 and homeodomain families. rCLAMPS learns a position-dependent amino-acid-nucleotide interaction energy model. To incorporate this model into the binding free energy calculation, we averaged the energy contributions across all occurrences of each amino-acid-nucleotide pair, which resulted in a 20-by-4 residue-type-specific energy matrix. This matrix is structurally analogous to the IDEA-trained energy model and can be directly integrated into the binding free energy calculations. As shown in Figure S6, Table S1, and Table S2, the IDEA model generally outperforms DBD-Hunter and rCLAMPS, demonstrating that it can achieve better predictive accuracy than both generic and family-specific knowledge-based models.

We also performed 10-fold cross-validation using the same HT-SELEX datasets, following the protocol described in the Methods Section Enhanced Modeling Prediction with SELEX Data. For each protein, we divided the entire dataset into 10 equal, randomly assigned folds. In each iteration, we used randomly selected 9 of the 10 folds as the training dataset and the remaining fold as the testing dataset. This process was repeated 10 times so that each fold served as the test set once. We then reported the average R² scores across these iterations to evaluate IDEA’s predictive performance. Our results are compared with the 1mer and 1mer+shape methods from,²⁶ the latest regression model that considers the shape of DNA with different sequences (Figure S7). This comparative analysis shows IDEA achieved higher predictive accuracy than the state-of-the-art sequence-based protein-DNA binding predictors for protein-DNA complexes that have available experimentally resolved structures.

Overall, these results demonstrate that IDEA can be used to predict the protein-DNA pairs in the absence of known binding sequence data, thus filling an important gap in protein-DNA predictions when experimental binding sequence data are unavailable.

Analysis of examples where IDEA fails to recognize strong DNA binders

We examine IDEA’s capability in identifying strong binders from the HT-SELEX dataset across 12 protein families.²⁶ The model successfully predicts 18 out of 22 protein-DNA systems, but the performance is reduced in 4 cases. Closer investigations revealed the source of these limitations. In some instances, only the protein from a different organism is available. For example, the PDX1 HT-SELEX data utilized the human PDX1 protein, but no human PDX1–DNA complex structure is available. Therefore, the mouse PDX1–DNA complex structure (PDB ID: 2H1K) was used for model training. Differences between model organisms may reduce predictive accuracy. A similar limitation applies to POU3F1, where an available mouse complex (PDB ID: 4Y60) was used to predict human protein–DNA interactions. Notably, DeepBind,²⁴ a sequence-based prediction tool, also failed to distinguish strong from weak binders when using the mouse POU3F1 protein (AUC score: 0.457), but this was corrected with the human POU3F1 protein (AUC score: 0.956).

IDEA also fails to fully resolve the binding preference of USF1. A closer examination of the HT-SELEX data reveals a lack of distinction among the sequences, as most sequences, including those with the lowest M-word (binding affinity) scores, contain the DNA-binding E-box sequence CACGTG. Therefore, USF1 represents a challenging example where the experimental data only consists of strong binders with limited variations in binding affinity, which likely results from differences in flanking sequences of the E-box motif.

Egr1 stands as a peculiar example. Whereas IDEA does not effectively distinguish between the strong and weak binders in the current HT-SELEX dataset, its predictions are consistent with other experimental datasets, including binding affinities measured by k-MITOMI⁴³ (Figure S8A, B), preferred binding sequences from protein-binding microarray, an earlier HT-SELEX experiment, and bacterial one-hybrid data.⁷ Therefore, further investigation of the current HT-SELEX data is needed to reconcile these differences.