A shape-constrained regression and wild bootstrap framework for reproducible drug synergy testing

High-throughput drug combination screens motivate computational methods to identify synergistic pairs, yet synergy is typically quantified by heuristic scores (Bliss, HSA, Loewe, ZIP) that provide no statistical inference and can be unstable or undefined when parametric dose-response fits fail. We present a nonparametric, assumption-light framework that defines interaction as the deviation from a monotone-additive null within a shared monotone model class. We fit a monotone surface by two-dimensional isotonic regression and a monotone-additive surface, compute an interaction surface, and summarize global interaction by a stable "interaction energy" statistic. A degrees-of-freedom-corrected wild bootstrap yields calibrated p-values for testing interaction in each dose-response matrix, enabling principled hit calling and multiple-testing control. On DrugCombDB, our method yields higher replicate concordance of interaction surfaces (median correlation 0.91 across 1,839 replicate pairs) than Bliss, HSA, Loewe, or ZIP (0.53-0.74), while avoiding the 20.9% Loewe and 3.6% ZIP failure rates. Because the fitted surface is generative, the method also predicts missing wells (median holdout RMSE 0.040 in viability units). By turning synergy scoring into statistically grounded outcomes (effect sizes with uncertainty), the framework provides more reliable targets for downstream machine learning models of combination response.