EPBDSurr

DNA breathing — the spontaneous, transient local opening of the double helix — modulates how transcription factors reach their binding sites, and the extended Peyrard-Bishop-Dauxois (EPBD) model captures it well. The problem is cost: pyDNA-EPBD runs Markov chain Monte Carlo and Langevin dynamics per sequence, putting genome-wide calculation out of reach at roughly 4,395 days on a 100-node CPU cluster.

Manish Bhattarai's group in the Theoretical Division at Los Alamos National Laboratory, with collaborators at George Mason and Brown, had already shown in Nucleic Acids Research that feeding EPBD features into a transformer improves transcription factor binding prediction. That model, EPBDxDNABERT-2, fuses DNABERT-2 sequence embeddings with biophysical features through cross-attention and a 690-way multi-label head, and beat a sequence-only baseline on more than 660 of 690 in vivo ChIP-seq binding events, with AUROC gains up to 9.6%. Its bottleneck was the simulator, not the network.

This work removes that bottleneck with EPBDSurr, a conditional denoising diffusion probabilistic model trained once on sequence-to-EPBD-feature pairs from pyDNA-EPBD, which then generates breathing features for new sequences with nothing re-fit per sequence.

#Key Features

• Surrogate for a physics simulator: The diffusion model learns the sequence-to-biophysics mapping that MCMC and Langevin dynamics compute explicitly, turning a months-long genome-wide calculation into a batched GPU inference pass.
• Sequence-conditioned generation: A transformer sequence embedder produces a fixed-length context vector conditioning every step of the reverse diffusion process, tying the biophysical output to the specific nucleotide sequence rather than sampling unconditionally.
• Distributional validation, not just error: Fidelity is assessed with Kolmogorov-Smirnov two-sample tests per feature and base-pair index, checking that generated and simulated samples share a distribution rather than only that their means agree.
• Drop-in for the downstream predictor: Surrogate features are consumed by the unchanged EPBDxDNABERT-2 architecture, so the cost of substituting them is measurable against ground-truth features.
• Chromosome-level splits: Chromosomes 1–6 and 10–22 train, chromosome 7 validates, and chromosomes 8 and 9 are held out, so generalization is tested across genomic regions rather than random windows.

#Technical Details

The surrogate has two components. The sequence embedder maps each nucleotide from the five-symbol alphabet to a learnable 256-dimensional vector, adds sinusoidal positional encodings, and passes the result through six transformer encoder layers with eight attention heads before attention pooling collapses it to one context vector. The denoiser is a UNet with four downsampling and four upsampling blocks at kernel size 3 and stride 2, with ReLU activations and dropout at 0.1; conditioning enters through adaptive group normalization modulated by the sequence embedding and a sinusoidal timestep embedding. Diffusion uses 1,000 timesteps under a linear noise schedule with beta rising from 1e-4 to 0.02. Training ran 100 epochs at batch size 128 with AdamW at an 8e-5 learning rate, cosine annealing, gradient clipping at norm 1.0 and a parameter EMA at decay 0.995, over more than a million 200-base-pair pairs.

Against simulated ground truth on 10,000 sequences the mean squared error is 0.0025, and KS tests fail to reject distributional equality for most feature and base-pair combinations at 95% confidence; the exception is base-flipping features at the 3.53 Å threshold, whose low signal magnitude the model does not capture. Substituting surrogate features into EPBDxDNABERT-2 without retraining gives AUROC 0.902 and AUPR 0.263; retraining recovers 0.940 and 0.293, against 0.949 and 0.326 on ground-truth features, and comparators DeepSEA (0.901/0.248), DanQ (0.931/0.295), TBiNet (0.947/0.333) and a fine-tuned DNABERT-2 (0.918/0.296). At batch size 8,192 the surrogate does in 8.3 minutes what pyDNA-EPBD takes 1,728 minutes to do; a 3-million-sequence ChIP-seq workload drops from 115.2 hours to 8.3 minutes, and the full genome from about 105,486 hours to 126.88 hours on an A100 cluster.

#Applications

The immediate use is genome-scale annotation: computing DNA breathing profiles across the human genome so that biophysically informed binding predictions run at sequence-only scale. That supports interpreting non-coding regulatory variants from genome-wide association studies, where the question is whether a substitution alters local helix stability enough to change occupancy. The paper demonstrates the necessary resolution on the adeno-associated virus P5 promoter, reproducing the coordinate-distance and base-flipping profiles that distinguish wild-type from a two-nucleotide mutant.

#Impact

The contribution is the demonstration that a narrow, well-posed biophysical simulation can be replaced by a conditional generative model at a fidelity high enough for the downstream task to barely notice — roughly one AUROC point against directly simulated features, for orders of magnitude less compute — a pattern that generalizes beyond DNA breathing. The surrogate systematically misses the weakest flipping features, the ones most tied to local melting and non-canonical structure, so the approximation is not uniform. The base EPBD-BERT code is BSD-3 licensed with a weights download script, but the surrogate's own code and weights are not in that repository, and the preprint gives no separate availability link.