Paper deep dive
Protein Counterfactuals via Diffusion-Guided Latent Optimization
Weronika Kłos, Sidney Bender, Lukas Kades
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 96%
Last extracted: 3/22/2026, 6:18:15 AM
Summary
MCCOP (Manifold-Constrained Counterfactual Optimization for Proteins) is a framework that generates minimal, biologically plausible sequence edits to flip protein property predictions. It operates in a continuous joint sequence-structure latent space, using a pretrained diffusion model as a manifold prior to ensure structural validity and sparsity, outperforming discrete baselines in GFP fluorescence, thermodynamic stability, and E3 ligase activity tasks.
Entities (6)
Relation Signals (4)
MCCOP → evaluatedon → GFP
confidence 98% · We evaluate MCCOP on three protein engineering tasks - GFP fluorescence rescue
MCCOP → evaluatedon → Ube4b
confidence 98% · We evaluate MCCOP on three protein engineering tasks - ... E3 ligase activity recovery
MCCOP → uses → DiMA
confidence 95% · We regularize the trajectory using DiMA (Meshchaninov et al., 2024) as an implicit manifold prior.
MCCOP → uses → CHEAP
confidence 95% · We map sequences to continuous representations using CHEAP (Lu et al., 2025).
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Deep learning models can predict protein properties with unprecedented accuracy but rarely offer mechanistic insight or actionable guidance for engineering improved variants. When a model flags an antibody as unstable, the protein engineer is left without recourse: which mutations would rescue stability while preserving function? We introduce Manifold-Constrained Counterfactual Optimization for Proteins (MCCOP), a framework that computes minimal, biologically plausible sequence edits that flip a model's prediction to a desired target state. MCCOP operates in a continuous joint sequence-structure latent space and employs a pretrained diffusion model as a manifold prior, balancing three objectives: validity (achieving the target property), proximity (minimizing mutations), and plausibility (producing foldable proteins). We evaluate MCCOP on three protein engineering tasks - GFP fluorescence rescue, thermodynamic stability enhancement, and E3 ligase activity recovery - and show that it generates sparser, more plausible counterfactuals than both discrete and continuous baselines. The recovered mutations align with known biophysical mechanisms, including chromophore packing and hydrophobic core consolidation, establishing MCCOP as a tool for both model interpretation and hypothesis-driven protein design. Our code is publicly available at this http URL.
Tags
Links
- Source: https://arxiv.org/abs/2603.10811v1
- Canonical: https://arxiv.org/abs/2603.10811v1
Trouble viewing inline? Open PDF directly →
Full Text
47,266 characters extracted from source content.
Expand or collapse full text
Protein Counterfactuals via Diffusion-Guided Latent Optimization Weronika Kłos1,2 Sidney Bender1,2 Lukas Kades3 1Machine Learning Group, Technische Universität Berlin, Berlin, Germany 2Berlin Institute for the Foundations of Learning and Data (BIFOLD) 3BASF Digital Solutions GmbH, Ludwigshafen am Rhein, Germany w.klos,s.bender@tu-berlin.de lukas.kades@basf.com Abstract Deep learning models can predict protein properties with unprecedented accuracy but rarely offer mechanistic insight or actionable guidance for engineering improved variants. When a model flags an antibody as unstable, the protein engineer is left without recourse: which mutations would rescue stability while preserving function? We introduce Manifold-Constrained Counterfactual Optimization for Proteins (MCCOP), a framework that computes minimal, biologically plausible sequence edits that flip a model’s prediction to a desired target state. MCCOP operates in a continuous joint sequence–structure latent space and employs a pretrained diffusion model as a manifold prior, balancing three objectives: validity (achieving the target property), proximity (minimizing mutations), and plausibility (producing foldable proteins). We evaluate MCCOP on three protein engineering tasks – GFP fluorescence rescue, thermodynamic stability enhancement, and E3 ligase activity recovery – and show that it generates sparser, more plausible counterfactuals than both discrete and continuous baselines. The recovered mutations align with known biophysical mechanisms, including chromophore packing and hydrophobic core consolidation, establishing MCCOP as a tool for both model interpretation and hypothesis-driven protein design. Our code is publicly available at github.com/weroks/mccop. 1 Introduction Deep learning has transformed computational protein science. Structure prediction models achieve near-experimental accuracy (Jumper et al., 2021; Abramson et al., 2024), protein language models capture evolutionary grammar (Lin et al., 2023; Team and others, 2024), and generative frameworks design novel folds from scratch (Watson et al., 2023; Ingraham et al., 2023). Yet these models remain oracles rather than guides: when a predictor flags a candidate as “aggregation-prone”, the engineer receives no indication of which mutations would resolve the problem. This paper addresses the need for algorithmic recourse: given a protein P predicted to lack a desired property ytargety_target, what is the minimal modification such that the prediction changes? This maps directly to counterfactual explanations (Wachter et al., 2017). Applied to a model of uncertain quality, counterfactuals expose reliance on spurious correlations; applied to a robust model, they generate testable hypotheses for wet-lab validation. Translating counterfactual methods to proteins introduces two fundamental challenges. First, the manifold constraint: Unlike images, proteins are governed by strict epistatic constraints – a single core mutation can abolish folding while a compensatory mutation restores it. Naive gradient optimization produces adversarial or invalid examples that satisfy the predictor but correspond to unfoldable proteins. Second, discreteness and geometry: Proteins are discrete sequences whose function emerges from continuous 3D geometry. Gradient-based methods require continuous relaxation, while naively treating them as sequences ignores spatial relationships: a mutation can be compensatory for another one only if the residues are proximal in 3D space, a property not directly apparent from the 1D sequence. We address both challenges with MCCOP, a gradient-based framework operating in a continuous joint sequence–structure embedding space that uses a pretrained diffusion model as a manifold prior. Our contributions are: 1. Framework. MCCOP combines predictor-guided gradient descent with diffusion-based manifold projection and gradient-sensitivity masking to produce sparse, valid, and plausible protein counterfactuals, without task-specific retraining of the generative model. 2. Quantitative evaluation. On three benchmarks, MCCOP achieves near-perfect success rates with 3–5× fewer mutations than discrete baselines and near-zero adversarial rates. 3. Mechanistic interpretability. MCCOP rediscovers known functional motifs and in several cases exactly recovers ground-truth counterfactual sequences from held-out test data. An overview of our approach is depicted in Figure 1. Figure 1: Overview of MCCOP. A non-fluorescent GFP variant is mapped to a continuous joint sequence–structure latent space via a pretrained autoencoder. After smoothing the classification boundary, the counterfactual embedding is optimized by alternating between (1) a sparse gradient step maximizing target class probability and (2) manifold projection using a pretrained diffusion model (DiMA). The counterfactual shown is a rediscovered sample from the held-out test set (red sticks denote the mutated residue). 2 Related Work Protein language models and embeddings. Models such as ESM-2 (Lin et al., 2023) and ESM-C (Team and others, 2024) learn unsupervised representations from millions of sequences. Recent multimodal embeddings go further: CHEAP (Lu et al., 2025) compresses ESMFold (Lin et al., 2023) activations into a joint sequence–structure representation whose decoder maps back to both amino acid sequences and atomistic coordinates. This bidirectional mapping is central to our approach. Diffusion and generative protein design. EvoDiff (Alamdari et al., 2023) and DiMA (Meshchaninov et al., 2024) apply diffusion to discrete sequences or continuous embeddings; RFdiffusion (Watson et al., 2023) and folding diffusion (Wu et al., 2024) operate in SE(3)SE(3) space (see Wang et al. (2025) for an overview). Most generative methods focus on conditional or unconditional sampling. MCCOP differs by using diffusion not for generation but as a regularizer within an optimization loop – conceptually an inversion of classifier guidance. Explainability in the protein domain. Prior work relies on attention visualization (Vig et al., 2020), feature attribution (Sibli et al., 2025; Dickinson and Meyer, 2022), gradient-based structure perturbation (Tan and Zhang, 2023), or sparse autoencoders applied to pLMs (Gujral et al., 2025). Unlike passive attribution, MCCOP provides active recourse: not just why a protein is predicted to fail, but how to rescue it. Counterfactual explanations. Counterfactuals, formalized for ML by Wachter et al. (2017), seek minimal input modifications that change a model’s output – a concept not to be confused with causal counterfactual inference in the structural causal model (SCM) sense (Pearl, 2009). Methods for tabular data (Mothilal et al., 2020; Russell, 2019) are well established. For high-dimensional inputs, diffusion-based approaches (DVCE (Augustin et al., 2022), DIME (Jeanneret et al., 2022), Diff-ICE (Pegios et al., 2025), FastDiME (Weng et al., 2024), ACE (Jeanneret et al., 2023)) generate on-manifold counterfactuals via guided denoising, with extensions to diverse sets (Bender et al., 2025; Bender and Morik, 2026), graphs (Bechtoldt and Bender, 2026; Chen et al., 2023), and text (Sarkar, 2024). GAN/VAE-based predecessors include DiVE (Rodriguez et al., 2021) and Diffeomorphic Counterfactuals (Dombrowski et al., 2024). To our knowledge, no prior work applies diffusion-guided counterfactual optimization to proteins. The closest biological relatives – latent fitness optimization (Ngo et al., 2024; Castro et al., 2022) – seek global optima rather than minimal edits and train task-specific generative models. 3 Methods We now describe each component of MCCOP: the latent representation, predictor smoothing, and the counterfactual optimization loop itself. 3.1 Problem Formulation Let ℳ⊂ℝL′×DM ^L × D denote the manifold of biologically plausible protein embeddings. Given a predictor fθ:ℳ→f_θ:M and an input embedding z0∈ℳz_0 with prediction y0=fθ(z0)y_0=f_θ(z_0), we seek: z∗=argminz∈ℳ[ℒtask(fθ(z),ytarget)+λd(z,z0)],z^*= _z [L_task(f_θ(z),\,y_target)+λ\,d(z,\,z_0) ], (1) where d enforces proximity to z0z_0 and z∈ℳz ensures plausibility. Without the manifold constraint, optimization yields adversarial examples (Dombrowski et al., 2024). We enforce it implicitly using the score function of a diffusion model trained on protein embeddings, whose denoising step acts as a projection Πℳ _M, interleaved with gradient steps on Eq. 1. 3.2 Latent Representation We map sequences S∈LS ^L to continuous representations z∈ℝL′×Dz ^L × D using CHEAP (Lu et al., 2025). The encoder ℰE compresses ESMFold activations into embeddings jointly capturing evolutionary and structural information. The decoder D maps z back to both a sequence S^=seq(z) S=D_seq(z) and backbone coordinates Ω^=struct(z) =D_struct(z), with near-perfect round-trip reconstruction (>>99% residue accuracy). Crucially, D is a position-wise MLP – each token S^i S_i depends only on ziz_i – enabling sequence-level sparsity via row-wise latent masking (§3.4.2). Both encoder and decoder are frozen throughout. 3.3 Predictor Smoothing Our framework is model-agnostic: any differentiable predictor on CHEAP embeddings can be used. As a test-bed we train a shallow MLP fθf_θ on flattened embeddings (architecture details in Appendix A). A non-smooth fθf_θ produces high-frequency gradients guiding optimization toward adversarial perturbations. Motivated by the observation of Bender et al. (2025), that smooth classifiers yield more reliable counterfactual optimization, we smooth fθf_θ via four complementary mechanisms: (1) spectral normalization (Miyato et al., 2018) on all linear layers; (2) Jacobian regularization (Jakubovitz and Giryes, 2018), penalizing ‖∇zfθ(z)‖F2\| _zf_θ(z)\|_F^2; (3) Softplus activations (β=1β=1); and (4) embedding-space adversarial augmentation via FGSM (Goodfellow et al., 2014), where perturbations decoding to the original sequence are added with the original label, teaching invariance to semantically null perturbations. As shown in Table 1, this reduces gradient norms by up to 4× while maintaining or improving AUROC. 3.4 Counterfactual Optimization Given embedding zorigz_orig with predicted class yorigy_orig, we seek z∗z^* such that fθ(z∗)=ytarget≠yorigf_θ(z^*)=y_target≠ y_orig, minimizing decoded mutations while staying on ℳM. Algorithm 1 summarizes the procedure. 3.4.1 Objective Function At step t, we minimize: ℒCF(zt)=log(1+exp(m−y~⋅fθ(zt)))⏟ℒmargin+λdist‖zt−zorig‖22⏟ℒproxL_CF(z_t)= (1+ (m- y· f_θ(z_t)) )_L_margin+ _dist \|z_t-z_orig\|_2^2_L_prox (2) where y~∈−1,+1 y∈\-1,+1\ is the signed target label, m>0m>0 is a confidence margin, and λdist _dist controls the proximity-validity trade-off. 3.4.2 Gradient-Based Sparsity Masking We compute per-position sensitivity si=‖∇ziℒCF‖2s_i=\| _z_iL_CF\|_2 and construct a binary mask selecting the top-k positions: Mi=[si≥s(k)].M_i=1 [s_i≥ s_(k) ]. (3) Gradients are applied only at masked positions; non-masked positions are hard-reset to zorigz_orig. Because D is position-wise, row-wise masking in latent space directly enforces sequence-space sparsity. The mask can alternatively be user-defined for constrained editing (e.g., fixing catalytic residues). 3.4.3 Manifold Projection We regularize the trajectory using DiMA (Meshchaninov et al., 2024) as an implicit manifold prior. At each step, we partially diffuse to noise level tdifft_diff, denoise to obtain Πϕ(zt′) _φ(z _t), and blend: zt+1=(1−α)zt′+αΠϕ(zt′),z_t+1=(1-α)\,z _t+α\, _φ(z _t), (4) where α∈[0,1]α∈[0,1] controls projection strength (α=0α=0: unconstrained; α=1α=1: full projection, which destabilizes optimization). We use α=0.3α=0.3 in practice (ablation in Appendix D). Algorithm 1 MCCOP: Manifold-Constrained Counterfactual Optimization for Proteins 0: Embedding zorigz_orig, predictor fθf_θ, diffusion projector Πϕ _φ, target label y~ y, sparsity k, projection strength α, margin m, learning rate η, max steps TmaxT_ , confidence threshold τ 1: z0←zorigz_0← z_orig 2: for t=0,1,…,Tmax−1t=0,1,…,T_ -1 do 3: Compute ℒCF(zt)L_CF(z_t) via Eq. 2 4: Compute per-position sensitivity: si=‖∇ziℒCF‖2s_i=\| _z_iL_CF\|_2 5: Construct top-k mask: Mi=[si≥s(k)]M_i=1[s_i≥ s_(k)] 6: Gradient step: zt′=zt−η⋅(M⊙∇zℒCF)z _t=z_t-η·(M _zL_CF) 7: Hard reset: zt′[i]←zorig[i]z _t[i]← z_orig[i] for all i where Mi=0M_i=0 8: Manifold projection: zt+1=(1−α)zt′+αΠϕ(zt′)z_t+1=(1-α)\,z _t+α\, _φ(z _t) 9: if σ(y~⋅fθ(zt+1))≥τσ( y· f_θ(z_t+1))≥τ and seq(zt+1)≠SorigD_seq(z_t+1)≠ S_orig then 10: return zt+1z_t+1 Early stopping: valid counterfactual found 11: end if 12: end for 13: return zTmaxz_T_ Return best attempt 3.5 Experimental Setup 3.5.1 Datasets We evaluate on three datasets with diverse physical origins (statistics in Appendix G): (1) TAPE Fluorescence (Sarkisyan et al., 2016; Rao et al., 2019): GFP homologs with bimodal fluorescence, binarized into bright/dark classes (optimize dark→ ). (2) TAPE Stability (Rocklin et al., 2017): proteolysis-based stability measurements; we remove the middle 33% quantile to create stable/unstable classes (optimize unstable→ ). (3) Ube4b Activity (Starita et al., 2013): ∼ 100k mutations in the U-box domain mapped to auto-ubiquitination activity; middle 33% removed, active/inactive classes defined by top/bottom quantiles (optimize inactive→ ). 3.5.2 Baselines We compare against: (1) Stochastic Hill Climbing: greedy random single-site mutations; (2) Genetic Algorithm: population-based evolution with edit-distance-penalized fitness; (3) Gradient Descent: unconstrained latent optimization without smoothing or manifold projection. Details in Appendix E. 3.5.3 Evaluation Metrics We assess validity and sparsity via success rate (fraction achieving target class), Hamming distance (number of mutations), and adversarial rate (fraction of successful counterfactuals corresponding to the same sequence). Structural plausibility is evaluated using ESM3-predicted pLDDT confidence and radius of gyration (RgR_g). Physicochemical plausibility is monitored via GRAVY hydrophobicity, instability index, and a binary solubility proxy. 4 Results We evaluate on complete test sets, excluding samples misclassified by the predictor. This results in n=2093n=2093, 22092209, and 26002600 samples for the stability, fluorescence, and activity datasets respectively. Results are mean ± std over three seeds. 4.1 Predictor Smoothing Improves Robustness Without Sacrificing Accuracy Table 1: Predictor AUROC and average L2L_2 gradient norm before and after smoothing (mean ± std, 3 seeds). Dataset AUROC (↑ ) Avg. L2L_2 norm (↓ ) Before After Before After Fluorescence 0.99±0.000.99± 0.00 0.99±0.000.99± 0.00 2.21±0.192.21± 0.19 1.10±0.101.10± 0.10 Stability 0.94±0.000.94± 0.00 0.98±0.010.98± 0.01 1.38±0.131.38± 0.13 0.36±0.080.36± 0.08 Activity 0.82±0.000.82± 0.00 0.93±0.010.93± 0.01 0.36±0.060.36± 0.06 0.33±0.040.33± 0.04 Table 1 shows that smoothing reduces gradient norms by up to 4× while maintaining or improving AUROC. The largest gain is on the activity dataset (AUROC: 0.82→ 0.93), likely because Jacobian regularization and adversarial augmentation reduce overfitting to noisy labels. 4.2 MCCOP Produces Valid and Sparse Counterfactuals Table 2: Success rate, adversarial rate, and edit distance (mean ± std, 3 seeds). Edit distance computed on successful counterfactuals only (confidence ≥0.95≥ 0.95, edit distance ≥1≥ 1). Discrete methods cannot produce adversarial examples by construction, so no values are bolded in this column. †Gradient Descent achieves 100% adversarial rate; edit distance is undefined. Dataset Method Success Rate (↑ ) Adv. Rate (↓ ) Edit Dist. (↓ ) Stability Genetic Algorithm 0.55±0.010.55± 0.01 0.00±0.000.00± 0.00 7.76±0.067.76± 0.06 Gradient Descent† 1.00±0.001.00± 0.00 1.00±0.001.00± 0.00 −- Stochastic Hill Climb 0.23±0.000.23± 0.00 0.00±0.000.00± 0.00 9.46±0.189.46± 0.18 MCCOP (ours) 1.00±0.00 17.77783pt6.44444pt1.00± 0.00 0.03±0.000.03± 0.00 2.32±0.01 17.77783pt6.44444pt2.32± 0.01 Fluorescence Genetic Algorithm 0.36±0.30 17.77783pt6.44444pt0.36± 0.30 0.00±0.000.00± 0.00 5.37±3.095.37± 3.09 Gradient Descent† 1.00±0.001.00± 0.00 1.00±0.001.00± 0.00 −- Stochastic Hill Climb 0.13±0.000.13± 0.00 0.00±0.000.00± 0.00 7.79±0.307.79± 0.30 MCCOP (ours) 0.19±0.000.19± 0.00 0.01±0.000.01± 0.00 1.37±0.01 17.77783pt6.44444pt1.37± 0.01 Activity Genetic Algorithm 0.17±0.150.17± 0.15 0.00±0.000.00± 0.00 6.24±3.676.24± 3.67 Gradient Descent† 1.00±0.001.00± 0.00 1.00±0.001.00± 0.00 −- Stochastic Hill Climb 0.03±0.000.03± 0.00 0.00±0.000.00± 0.00 10.91±0.0210.91± 0.02 MCCOP (ours) 1.00±0.00 17.77783pt6.44444pt1.00± 0.00 0.02±0.020.02± 0.02 2.46±0.33 17.77783pt6.44444pt2.46± 0.33 Table 2 reveals three key findings. (1) Unconstrained gradient optimization is entirely adversarial: every counterfactual decodes to the original sequence, confirming exploitable high-frequency artifacts and validating our smoothing and projection pipeline. (2) MCCOP achieves high success with minimal edits: 100% success on stability and activity with 2.3–2.5 mutations versus 6.2–10.9 for discrete baselines. MCCOP reaches early stopping after a median of 2–10 steps, while hill climbing exhausts the budget in >>95% of cases. (3) Fluorescence is harder: MCCOP’s 19% success rate reflects the requirement for precise chromophore geometry potentially exceeding our sparsity budget (k=5k=5), yet successful counterfactuals are the sparsest (1.4 mutations) with near-zero adversarial rate. Edit distances for MCCOP and the genetic algorithm are tunable via k/λdist _dist and fitness weighting, respectively; MCCOP’s advantage lies in the favorable trade-off between success rate, sparsity, and plausibility. 4.3 Structural and Physicochemical Plausibility Figure 2: Physicochemical plausibility across benchmarks (columns: pLDDT, GRAVY, instability index, RgR_g; rows: fluorescence, activity, stability). MCCOP (orange) closely matches the original distribution (gray); discrete baselines show broader shifts. Statistical comparisons via Kruskal-Wallis/Dunn’s tests with Benjamini-Hochberg correction: MCCOP achieves significantly higher pLDDT than both baselines across all tasks (adjusted p<0.02p<0.02). Figure 2 shows that MCCOP counterfactuals are nearly indistinguishable from the original distribution across all metrics, occasionally shifting toward more favorable values. Discrete baselines introduce broader shifts, especially in hydrophobicity and instability index, as they explore sequence space without structural priors. A controlled comparison at fixed edit distance (Appendix C) confirms these trends. 4.4 MCCOP Rediscovers Known Biophysical Mechanisms Figure 3: Per-residue mutation frequency for fluorescence (A) and Ube4b activity (B). MCCOP (blue) concentrates mutations in functionally relevant regions – chromophore-proximal residues for GFP, E2-binding interface for Ube4b – while baselines distribute mutations nearly uniformly. Shaded regions: known functional motifs (Sarkisyan et al., 2016; Starita et al., 2013). GFP fluorescence. MCCOP concentrates mutations in the chromophore-proximal region (residues 63–69) and β-barrel strands forming the chromophore cavity (Figure 3A), consistent with the requirement for tight packing to suppress non-radiative decay (Sarkisyan et al., 2016). A small number of distal mutations (e.g., residues 181, 216) may represent novel compensatory interactions or predictor artifacts, requiring experimental follow-up. Ube4b activity. Mutations cluster at the E2-binding interface (residues 66–71; Figure 3B), through which Ube4b recruits UbcH5c for ubiquitin transfer (Starita et al., 2013). Thermodynamic stability. The stability dataset spans diverse topologies (Rocklin et al., 2017), so no universal residue positions dominate. However, MCCOP frequently targets core-facing residues, suggesting hydrophobic core consolidation as a general stabilization strategy (Appendix F). Recovery of ground-truth counterfactuals. MCCOP exactly recovers existing opposite-label sequences in 16 (fluorescence), 18 (activity), and 4 (stability) cases – several from the held-out test set. Figure 4 shows structural alignments confirming that recovered mutations localize to functionally relevant regions. Figure 4: Structural alignments between original (gray) and counterfactual (colored) proteins for rediscovered ground-truth examples. (A) GFP: mutations near the chromophore. (B) Stability: core-facing mutations. (C) Ube4b: E2-binding interface mutations. Structures predicted by ESM3. 5 Discussion MCCOP generates sparse, on-manifold counterfactual explanations achieving near-perfect success rates with 1–3 mutations on average (vs. 8–11 for discrete baselines) while maintaining structural and physicochemical plausibility. Recovered mutations align with established biophysical mechanisms, suggesting that the underlying predictors have learned meaningful sequence-function relationships. Explanation versus editing. MCCOP’s primary goal is model interpretation, not direct engineering. A counterfactual is only as trustworthy as the predictor it explains: if the model has learned spurious correlations, the counterfactual reflects them faithfully – which is itself diagnostic. When the predictor is robust, MCCOP’s outputs become candidates for experimental validation. From correlation to causation. Our framework identifies correlational, not causal, relationships. Establishing true causal links requires interventional experiments, but MCCOP’s sparse suggestions (2 mutations vs. thousands of directed-evolution variants) are directly amenable to such follow-up. Limitations. (1) Plausibility evaluation relies on computational proxies (ESM3 pLDDT, RgR_g, physicochemical indices) rather than experimental validation. (2) The CHEAP encoder–decoder introduces reconstruction error that may produce artifacts for proteins distant from ESMFold’s training distribution. (3) We evaluate only binary tasks; extending to continuous regression targets requires replacing the margin loss with MSE or quantile losses. On the manifold and smoothness assumptions. Two assumptions embedded in our framework deserve scrutiny. First, MCCOP operates in a continuous latent space under the implicit premise that plausible protein sequences concentrate near a low-dimensional manifold. The same manifold hypothesis is routinely invoked in computer vision, where natural images are assumed to populate a thin subspace of pixel space, yet to the best of our knowledge no formal proof of this assertion exists for images or for proteins. Fefferman et al. develop statistical tests for the hypothesis but do not establish it for any natural data distribution (Fefferman et al., 2016); empirically, the evidence is consistent with data concentrating on disconnected clusters or “blobs” rather than a single smooth manifold (Bengio et al., 2013). For proteins, the situation is arguably more fraught: functional sequences are constrained by folding, stability, and epistasis, producing a viable sequence space that may be fragmented and topologically complex rather than smoothly connected. Second, MCCOP’s Gaussian smoothing of the latent space presupposes that the underlying sequence–function mapping varies smoothly, so that local perturbations yield gradual changes in the predicted phenotype. However, protein fitness landscapes are known to be rugged: higher-order epistasis creates abrupt fitness transitions even between sequences that differ by a single residue (Weinreich et al., 2006; Sarkisyan et al., 2016), and there is no a priori reason to expect the predictor’s decision surface, which reflects these landscapes, to be smooth either. An alternative to smoothing might be signal filtering tuned to a desired frequency, which would suppress high-frequency noise without globally flattening the landscape; we opted for Gaussian smoothing as a pragmatic engineering choice that made the gradient-based optimization tractable, rather than as a theoretically motivated operation. We flag these points not because they invalidate the results – MCCOP’s strong empirical performance suggests the assumptions are serviceable in practice – but because they circumscribe the regime in which the method’s outputs should be trusted and highlight opportunities for more principled geometric and spectral approaches in future work. Future directions. (1) Multi-objective counterfactuals: jointly optimizing stability and binding affinity by combining predictor gradients. (2) Experimental validation: synthesizing top-ranked variants for closed-loop validation. (3) Diverse counterfactual sets (Mothilal et al., 2020; Bender and Morik, 2026): revealing alternative mutational strategies and enriching fitness landscape understanding. Acknowledgments We would like to thank Marvin Sextro for many valuable pieces of advice and proof-reading, as well as Klaus-Robert Müller, Adrian Hill, and Stefan Chmiela for interesting and fruitful discussions. We also would like to thank Dominik Kühne for maintaining our HPC cluster hydra and being always there to help in case of technical difficulties. We used GitHub Copilot for assistance with code development and editing of paper text. All AI-generated content was reviewed, verified, and revised by the authors, who take full responsibility for the final manuscript. This work was supported by the German Ministry for Education and Research (BMBF) under Grant 01IS18037A, and by BASLEARN – TU Berlin/BASF Joint Laboratory, co-financed by TU Berlin and BASF SE. References J. Abramson, J. Adler, J. Dunger, R. Evans, T. Green, A. Pritzel, O. Ronneberger, L. Willmore, A. J. Ballard, J. Bambrick, et al. (2024) Accurate structure prediction of biomolecular interactions with alphafold 3. Nature 630 (8016), p. 493–500. Cited by: §1. S. Alamdari, N. Thakkar, R. Van Den Berg, N. Tenenholtz, R. Strome, A. M. Moses, A. X. Lu, N. Fusi, A. P. Amini, and K. K. Yang (2023) Protein generation with evolutionary diffusion: sequence is all you need. BioRxiv, p. 2023–09. Cited by: §2. M. Augustin, V. Boreiko, F. Croce, and M. Hein (2022) Diffusion visual counterfactual explanations. Advances in Neural Information Processing Systems 35, p. 364–377. Cited by: §2. D. Bechtoldt and S. Bender (2026) Graph diffusion counterfactual explanation. ESANN. Cited by: §2. S. Bender, J. Herrmann, K. Müller, and G. Montavon (2025) Towards desiderata-driven design of visual counterfactual explainers. Pattern Recognition. Cited by: §2, §3.3. S. Bender and M. Morik (2026) Visual disentangled diffusion autoencoders. arXiv preprint arXiv:2601.21851. Cited by: §2, §5. Y. Bengio, A. Courville, and P. Vincent (2013) Representation learning: a review and new perspectives. IEEE transactions on pattern analysis and machine intelligence 35 (8), p. 1798–1828. Cited by: §5. E. Castro, A. Godavarthi, J. Rubinfien, K. B. Givechian, D. Bhaskar, and S. Krishnaswamy (2022) ReLSO: a transformer-based model for latent space optimization and generation of proteins. arXiv preprint arXiv:2201.09948. Cited by: §2. J. Chen, S. Wu, A. Gupta, and R. Ying (2023) D4explainer: in-distribution explanations of graph neural network via discrete denoising diffusion. Advances in Neural Information Processing Systems 36, p. 78964–78986. Cited by: §2. Q. Dickinson and J. G. Meyer (2022) Positional shap (poshap) for interpretation of machine learning models trained from biological sequences. PLOS Computational Biology 18 (1), p. e1009736. Cited by: §2. A. Dombrowski, J. E. Gerken, K. Müller, and P. Kessel (2024) Diffeomorphic counterfactuals with generative models. IEEE Transactions on Pattern Recognition and Machine Intelligence 46 (5), p. 3257–3274. Cited by: §2, §3.1. C. Fefferman, S. Mitter, and H. Narayanan (2016) Testing the manifold hypothesis. Journal of the American Mathematical Society 29 (4), p. 983–1049. Cited by: §5. I. J. Goodfellow, J. Shlens, and C. Szegedy (2014) Explaining and harnessing adversarial examples. arXiv preprint arXiv:1412.6572. Cited by: §3.3. O. Gujral, M. Bafna, E. Alm, and B. Berger (2025) Sparse autoencoders uncover biologically interpretable features in protein language model representations. Proceedings of the National Academy of Sciences 122 (34), p. e2506316122. Cited by: §2. J. B. Ingraham, M. Baranov, Z. Costello, K. W. Barber, W. Wang, A. Ismail, V. Frappier, D. M. Lord, C. Ng-Thow-Hing, E. R. Van Vlack, et al. (2023) Illuminating protein space with a programmable generative model. Nature 623 (7989), p. 1070–1078. Cited by: §1. D. Jakubovitz and R. Giryes (2018) Improving dnn robustness to adversarial attacks using jacobian regularization. In Proceedings of the European conference on computer vision (ECCV), p. 514–529. Cited by: item 2, §3.3. G. Jeanneret, L. Simon, and F. Jurie (2022) Diffusion models for counterfactual explanations. In Proceedings of the Asian Conference on Computer Vision, p. 858–876. Cited by: §2. G. Jeanneret, L. Simon, and F. Jurie (2023) Adversarial counterfactual visual explanations. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 16425–16435. Cited by: §2. J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko, et al. (2021) Highly accurate protein structure prediction with alphafold. nature 596 (7873), p. 583–589. Cited by: §1. Z. Lin, H. Akin, R. Rao, B. Hie, Z. Zhu, W. Lu, N. Smetanin, R. Verkuil, O. Kabeli, Y. Shmueli, et al. (2023) Evolutionary-scale prediction of atomic-level protein structure with a language model. Science 379 (6637), p. 1123–1130. Cited by: §1, §2. A. X. Lu, W. Yan, K. K. Yang, V. Gligorijevic, K. Cho, P. Abbeel, R. Bonneau, and N. C. Frey (2025) Tokenized and continuous embedding compressions of protein sequence and structure. Patterns 6 (6). Cited by: §2, §3.2. V. Meshchaninov, P. Strashnov, A. Shevtsov, F. Nikolaev, N. Ivanisenko, O. Kardymon, and D. Vetrov (2024) Diffusion on language model encodings for protein sequence generation. arXiv preprint arXiv:2403.03726. Cited by: §2, §3.4.3. T. Miyato, T. Kataoka, M. Koyama, and Y. Yoshida (2018) Spectral normalization for generative adversarial networks. arXiv preprint arXiv:1802.05957. Cited by: item 1, Appendix A, §3.3. R. K. Mothilal, A. Sharma, and C. Tan (2020) Explaining machine learning classifiers through diverse counterfactual explanations. In Proceedings of the 2020 conference on fairness, accountability, and transparency, p. 607–617. Cited by: §2, §5. N. K. Ngo, T. V. Tran, V. T. Duy Nguyen, and T. S. Hy (2024) Latent-based directed evolution accelerated by gradient ascent for protein sequence design. bioRxiv, p. 2024–04. Cited by: §2. J. Pearl (2009) Causality. Cambridge university press. Cited by: §2. P. Pegios, M. Lin, N. Weng, M. B. S. Svendsen, Z. Bashir, S. Bigdeli, A. N. Christensen, M. Tolsgaard, and A. Feragen (2025) Diffusion-based iterative counterfactual explanations for fetal ultrasound image quality assessment. In International Workshop on Advances in Simplifying Medical Ultrasound, p. 174–184. Cited by: §2. R. Rao, N. Bhattacharya, N. Thomas, Y. Duan, P. Chen, J. Canny, P. Abbeel, and Y. Song (2019) Evaluating protein transfer learning with tape. Advances in neural information processing systems 32. Cited by: §3.5.1. G. J. Rocklin, T. M. Chidyausiku, I. Goreshnik, A. Ford, S. Houliston, A. Lemak, L. Carter, R. Ravichandran, V. K. Mulligan, A. Chevalier, et al. (2017) Global analysis of protein folding using massively parallel design, synthesis, and testing. Science 357 (6347), p. 168–175. Cited by: §3.5.1, §4.4. P. Rodriguez, M. Caccia, A. Lacoste, L. Zamparo, I. Laradji, L. Charlin, and D. Vazquez (2021) Beyond trivial counterfactual explanations with diverse valuable explanations. In Proceedings of the IEEE/CVF International Conference on Computer Vision, p. 1056–1065. Cited by: §2. C. Russell (2019) Efficient search for diverse coherent explanations. In Proceedings of the conference on fairness, accountability, and transparency, p. 20–28. Cited by: §2. A. Sarkar (2024) Large language models cannot explain themselves. arXiv preprint arXiv:2405.04382. Cited by: §2. K. S. Sarkisyan, D. A. Bolotin, M. V. Meer, D. R. Usmanova, A. S. Mishin, G. V. Sharonov, D. N. Ivankov, N. G. Bozhanova, M. S. Baranov, O. Soylemez, et al. (2016) Local fitness landscape of the green fluorescent protein. Nature 533 (7603), p. 397–401. Cited by: §3.5.1, Figure 3, §4.4, §5. S. A. Sibli, V. P. Panagiotou, and C. Makris (2025) Enhancing protein structure predictions: deepshap as a tool for understanding alphafold2. Expert Systems with Applications, p. 127853. Cited by: §2. L. M. Starita, J. N. Pruneda, R. S. Lo, D. M. Fowler, H. J. Kim, J. B. Hiatt, J. Shendure, P. S. Brzovic, S. Fields, and R. E. Klevit (2013) Activity-enhancing mutations in an e3 ubiquitin ligase identified by high-throughput mutagenesis. Proceedings of the National Academy of Sciences 110 (14), p. E1263–E1272. Cited by: §3.5.1, Figure 3, §4.4. J. Tan and Y. Zhang (2023) Explainablefold: understanding alphafold prediction with explainable ai. In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, p. 2166–2176. Cited by: §2. E. Team et al. (2024) ESM cambrian: revealing the mysteries of proteins with unsupervised learning. Evolutionary Scale Website https://w. evolutionaryscale. ai/blog/esm-cambrian. Cited by: §1, §2. J. Vig, A. Madani, L. R. Varshney, C. Xiong, R. Socher, and N. F. Rajani (2020) Bertology meets biology: interpreting attention in protein language models. arXiv preprint arXiv:2006.15222. Cited by: §2. S. Wachter, B. Mittelstadt, and C. Russell (2017) Counterfactual explanations without opening the black box: automated decisions and the gdpr. Harv. JL & Tech. 31, p. 841. Cited by: §1, §2. C. Wang, S. Alamdari, C. Domingo-Enrich, A. P. Amini, and K. K. Yang (2025) Toward deep learning sequence–structure co-generation for protein design. Current Opinion in Structural Biology 91, p. 103018. Cited by: §2. J. L. Watson, D. Juergens, N. R. Bennett, B. L. Trippe, J. Yim, H. E. Eisenach, W. Ahern, A. J. Borst, R. J. Ragotte, L. F. Milles, et al. (2023) De novo design of protein structure and function with rfdiffusion. Nature 620 (7976), p. 1089–1100. Cited by: §1, §2. D. M. Weinreich, N. F. Delaney, M. A. DePristo, and D. L. Hartl (2006) Darwinian evolution can follow only very few mutational paths to fitter proteins. science 312 (5770), p. 111–114. Cited by: §5. N. Weng, P. Pegios, E. Petersen, A. Feragen, and S. Bigdeli (2024) Fast diffusion-based counterfactuals for shortcut removal and generation. In European Conference on Computer Vision, p. 338–357. Cited by: §2. K. E. Wu, K. K. Yang, R. van den Berg, S. Alamdari, J. Y. Zou, A. X. Lu, and A. P. Amini (2024) Protein structure generation via folding diffusion. Nature communications 15 (1), p. 1059. Cited by: §2. Appendix A Predictor Training and Smoothing Details Architecture. The property predictor fθf_θ is a three-layer MLP with hidden dimensions [512, 256], each followed by spectral normalization (Miyato et al., 2018) and Softplus activation (β=1β=1). The final layer outputs a single logit. Input embeddings are flattened resulting in an input dimension of sequence length L times embedding dimension D (masking padding tokens) before being passed to the MLP. Training protocol. We train on 80% of each dataset, reserving 10% for validation and 10% for testing, stratified by label. We use Adam with a learning rate of 1×10−51× 10^-5 and a dropout rate of 0.30.3. Early stopping is applied with a patience of 5 epochs based on validation AUROC. Smoothing mechanisms. Further details concerning the smoothing mechanisms include: 1. Spectral normalization (Miyato et al., 2018): applied to all linear layers, constraining the Lipschitz constant of each layer to approximately 1. 2. Jacobian regularization (Jakubovitz and Giryes, 2018): we add a penalty term λJ‖∇zfθ(z)‖F2 _J\| _zf_θ(z)\|_F^2 to the training loss, with λJ=10−3 _J=10^-3. The Frobenius norm is estimated via a Hutchinson trace estimator with 5 random projections per batch for computational efficiency. 3. Adversarial data augmentation: for each training sample (zi,yi)(z_i,y_i), we generate an adversarial embedding ziadvz_i^adv via an FGSM attack (ϵ=0.01ε=0.01 in embedding space) targeting the opposite class. Adversarial samples that decode to the same amino acid sequence as the original (i.e., (ziadv)=(zi)D(z_i^adv)=D(z_i)) are added to the training set with the original label yiy_i, teaching the model to be invariant to semantically null perturbations. Smoothness quantification. We report the average L2L_2 gradient norm z∼test[‖∇zfθ(z)‖2]E_z _test[\| _zf_θ(z)\|_2] computed over the full test set. Lower values indicate a smoother decision boundary. Appendix B Computational Cost Analysis Figure 5 reports the average wall-clock time per sample for MCCOP and the two discrete baselines across all three benchmarks. The genetic algorithm is the most expensive method by roughly an order of magnitude due to its population-based evaluation. MCCOP and stochastic hill climbing have comparable per-sample execution times. An important caveat applies to this comparison. The discrete baselines operate in sequence space but must evaluate candidates using the same embedding-space predictor; every proposed mutation therefore requires a full re-encoding through the CHEAP encoder (backed by ESMFold), which constitutes 97% and 94% of total computation time for hill climbing and the genetic algorithm respectively. This overhead is not intrinsic to those algorithms but arises from the requirement of a shared evaluation protocol. MCCOP, by contrast, operates natively in embedding space and avoids this round-trip entirely. Its dominant cost is the diffusion-based manifold projection, which accounts for 99% of computation time. For performance-critical applications we would recommend executing the diffusion-based projection step only every n optimization steps. Figure 5: Average wall-clock time per sample across the complete test set. The genetic algorithm is the most expensive method. MCCOP and stochastic hill climbing have comparable execution times, though their computational profiles differ: discrete baselines spend >>94% of time on re-encoding candidate sequences, while MCCOP spends 99% on diffusion-based manifold projection. Appendix C Controlled Edit Distance Comparison To ensure a fair comparison across methods, we filter all successful counterfactuals to those with exactly three mutations (edit distance = 3), which represents the bin with the highest overlap across MCCOP, the genetic algorithm, and stochastic hill climbing. Figure 6 shows the same physicochemical property distributions as Figure 2 in the main text, restricted to this subset. The trends observed in the main text are preserved: MCCOP-generated counterfactuals remain within the distribution of the original test set for pLDDT, GRAVY, instability index, and radius of gyration, whereas the discrete baselines show broader deviations, particularly in instability index and GRAVY. Figure 6: Physicochemical property distributions for counterfactuals with exactly 3 mutations. MCCOP closely tracks the original test set distribution, while the genetic algorithm and stochastic hill climbing show greater deviation, particularly for GRAVY and instability index. Appendix D Hyperparameter Sensitivity Table 3 lists the primary hyperparameters for MCCOP and the values used. We use the same set of hyperparameters across all three datasets. Table 3: MCCOP hyperparameters and their values across benchmarks. Symbol Description Fluorescence Stability Activity k Top-k masked positions 5 5 5 λdist _dist L2L_2 distance weight 0.1 0.1 0.1 m Margin in ℒmarginL_margin 2.2 2.2 2.2 α Projection strength 0.3 0.3 0.3 tdifft_diff Diffusion noise level 100 100 100 η Learning rate 0.50.5 0.50.5 0.50.5 TmaxT_ Max optimization steps 50 50 50 We perform an ablation over all combinations of the previously mentioned smoother components (including no smoothing) as well as the manifold projection step and the masking value k. No smoothing, no masking and no manifold projection corresponds exactly to the gradient descent baseline. Lower k values proved to be too restrictive and caused low success rates while higher ones provided only a marginal increase in success rate. No smoothing significantly increases adversarial rates, while no manifold projection significantly reduces pLDDT scores of generated counterfactuals. Appendix E Baselines Implementation Details We compare our method against three baseline counterfactual explanation strategies, each operating over protein sequences and their corresponding embeddings. All baselines share a common interface and are evaluated using the same predictor and confidence threshold (τ=0.95τ=0.95 by default). Below we describe each baseline along with its hyperparameters. E.1 Gradient Descent This baseline performs standard gradient descent directly in the continuous embedding space. Given an input embedding x, a differentiable copy ′x is optimized to maximize the predictor’s probability of the target (flipped) class via binary cross-entropy loss. The Adam optimizer is used to update ′x over a fixed number of steps. At each step, the candidate counterfactual with the highest confidence toward the target class is retained. Table 4: Hyperparameters for the gradient descent baseline. Hyperparameter Symbol Value Learning rate η 1×10−21× 10^-2 Gradient steps T 5050 Confidence threshold τ 0.950.95 Optimizer – Adam Loss function ℒL BCEWithLogitsLoss Notably, this baseline operates entirely in embedding space and does not enforce any manifold constraints or discrete sequence validity, making it a purely continuous relaxation approach. E.2 Random Mutation This baseline performs a stochastic hill-climbing search in discrete sequence space. At each step, a single random point mutation is applied to each unsolved sequence: a uniformly random position is selected and replaced with a uniformly random amino acid from the standard 20-letter alphabet. The mutated sequence is re-encoded into embedding space using a lightweight encoder, and the predictor evaluates the new embedding. If the target-class confidence improves, the mutation is accepted; otherwise, the sequence reverts to the previous best. The process repeats for a fixed number of steps. Table 5: Hyperparameters for the Random Mutation baseline. Hyperparameter Symbol Value Maximum steps T 5050 Confidence threshold τ 0.950.95 Amino acid alphabet A Standard 20 Mutations per step – 11 E.3 Genetic Algorithm This baseline employs a population-based evolutionary strategy operating in discrete sequence space. For each input sequence, an initial population is constructed by applying random mutations to the original. At each generation, individuals are evaluated by encoding them into embedding space and computing a fitness score defined as the predictor’s target-class confidence, optionally penalized by the Hamming distance to the original sequence: f()=conf()−λ⋅dH(,orig),f(s)=conf(s)-λ· d_H(s,s_orig), (5) where conf()conf(s) is the predictor’s confidence on the target class for sequence s, dHd_H denotes the Hamming distance, and λ is the edit distance penalty. Selection uses tournament selection with tournament size 3. The top 20% of the population is preserved as elites. Offspring are generated via single-point crossover and random point mutation (1–2 mutations per offspring). Evolution proceeds for a fixed number of generations or until all samples exceed the confidence threshold. Table 6: Hyperparameters for the Genetic Algorithm baseline. Hyperparameter Symbol Value Population size N 4040 Generations G 3030 Crossover rate pcp_c 0.50.5 Edit distance penalty λ 0.020.02 Confidence threshold τ 0.950.95 Elite fraction – 20%20\% Tournament size k 33 Mutations per offspring – 11–22 Maximum batch size – 88 Appendix F Additional Structural Visualizations We provide two additional structure visualizations for the stability dataset as we could not verify hydrophobic core packing by investigating mutation frequencies per residue. Figure 7: Structural alignment of original (gray) and counterfactual (cyan) stability variants across three topologies. Core-facing mutated residues shown as sticks. Appendix G Dataset Statistics and Preprocessing Table 7: Dataset statistics after preprocessing. Dataset Sequences Positive Negative Avg. Length Binarization Fluorescence 54,025 30,697 23,328 236.96 Bimodal split Stability 45,901 22,694 23,207 45.06 Remove middle 33% Activity 60,692 30,293 30,399 102 Remove middle 33% For the fluorescence dataset, we exploit the natural bimodality of the log-fluorescence distribution and determine the optimal threshold using Otsu’s method. For the stability and activity datasets, removing the middle tercile ensures a clear margin between classes, reducing label noise near the decision boundary. All embeddings are computed using the CHEAP encoder with ESMFold as the backbone, producing representations of dimension D=1024D=1024 with no compression along the length dimension.