{"id":"625a679f-3ee2-45ca-9220-c26f0da9af8f","arxiv_id":"2607.22338","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":6,"one_line_summary":"Hierarchical GP inference with Monte-Carlo-sampled hyperparameters yields free-energy uncertainty estimates that track reconstruction error across data-ablation tests, unlike fixed-hyperparameter GP and umbrella integration baselines.","lead":"The paper builds a hierarchical Gaussian-process model that treats its own noise and length-scale settings as unknowns sampled by Monte Carlo, producing free-energy uncertainty estimates that change with the amount of simulation data. A reader doing umbrella sampling or metadynamics would care because these uncertainties could tell when a free-energy profile is actually converged, whereas the authors show conventional GP baselines give nearly constant or misleading error bars","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The headline claim that hierarchical GP uncertainty 'tracks' error is not yet quantitatively established: umbrella-sampling evidence is visual, and no pointwise calibration or coverage analysis is reported, so the acknowledged stationary-kernel assumption remains an untested threat to the UQ claim.","rationale":"The central empirical claim is comparative and about UQ usefulness; it would be true only if the hierarchical model's intervals are informative about error. The paper gives good evidence that SD changes with data and that MAP underestimates hyperparameter uncertainty, but it never measures the calibration of those intervals. The stationary SE kernel plus scalar noise is the obvious place where misspecification can enter: if the kernel is wrong, noise hyperparameters absorb bias and intervals can be locally misleading while average SD still correlates with RMSE. This is not an accusation of error; the authors explicitly flag the stationary-kernel limitation. My concern is that the validation tools used do not establish the central claim strongly enough. A pointwise calibration/coverage analysis on the existing 100-cell ablation grid is cheap because code/data are provided and would settle the issue. I agree with the reader's conditional verdict: the method is promising and well-engineered, but the headline claim should await calibration evidence or be restricted to 'SD responds to data amount' rather than 'uncertainty tracks error.'","tokens_in":21479,"tokens_out":12598,"duration_ms":120544,"concrete_test":"Using the provided R9 umbrella data, reproduce the 100-condition T×W ablation grid. For each condition and each CV bin, compute the standardized residual z_i = (A_WHAM(ξ_i) − μ_GP(ξ_i)) / σ_GP(ξ_i). Report the empirical coverage of the 68% and 95% posterior intervals, averaged over window subsets, and plot coverage versus CV position and data condition. Repeat for the fixed-GP, MAP, and UI baselines. If the hierarchical GP's 68%/95% coverage is close to nominal across all cells, the concern is resolved; if it undercovers in barrier regions or low-data cells, the stationary/homoskedastic prior is the likely cause and the 'tracking' claim should be restricted.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section 5 claims hierarchical GP inference was 'the only approach whose uncertainty estimates consistently adapted to the amount and quality' of the data. For this claim to hold, predictive intervals must not just shrink with data; their width must track the actual reconstruction error. The umbrella-sampling support (Fig. 3) is visual: averaged RMSE and averaged SD heatmaps, with no pointwise calibration, coverage, or standardized-residual analysis. The metadynamics support (Sec. 3.4) is a Pearson rho=0.95 computed on a small number of trajectory-length aggregates. Meanwhile the GP uses a stationary SE kernel (Eq. 23) and scalar inferred noise, a mismatch the authors explicitly acknowledge: 'there is no physical reason why a free energy profile should exhibit uniform correlation structure throughout an entire CV domain.' If true length scales vary near barriers/basins, the inferred noise and length scale can absorb the misspecification and produce intervals that are locally over- or under-confident even if the average SD correlates with RMSE. Because the evaluation never checks interval coverage, the central claim is conditional, not established.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper develops a hierarchical Bayesian Gaussian process (GP) framework for uncertainty quantification in free energy calculations. A zero-mean GP with squared-exponential kernel models the free energy profile from joint histogram and derivative observations, using a leave-one-out pseudo-likelihood and Hamiltonian Monte Carlo (HMC-NUTS) to sample the hyperposterior of kernel length scale, amplitude, and noise hyperparameters. Predictive uncertainty is obtained by propagating hyperparameter samples through the GP predictive equations, thereby including hyperparameter uncertainty. The method is applied to umbrella sampling data for an R9 peptide–membrane system and to extended Lagrangian metadynamics data for phenol–membrane interactions, with data-ablation studies comparing the proposed method against umbrella integration with block averaging, a fixed-hyperparameter GP, and a MAP-optimized GP. The central claim is that the hierarchical GP is the only approach whose uncertainty estimates consistently adapt to the amount and quality of simulation data, as assessed by visual heatmaps of average RMSE and average predictive standard deviation and by a Pearson correlation in the metadynamics case.","tokens_in":21838,"tokens_out":3992,"duration_ms":38346,"significance":"If the central claim is established, the method could provide a practically useful convergence signal for enhanced-sampling free energy calculations, with direct implications for automated workflows and for avoiding wasted simulation effort. The paper contains several strengths: the GP machinery is standard and correctly assembled; the HMC diagnostics (Section C) are thorough, including R-hat, ESS, and divergence counts across multiple data regimes; a hyperprior sensitivity analysis (Section B.2) is reported; and code and data are made available on GitHub and Zenodo. The authors also explicitly acknowledge the key limitation of using a stationary covariance kernel. However, the headline validation is currently incomplete: the umbrella-sampling evidence is qualitative, the metadynamics evidence rests on a single aggregate correlation, and no coverage or calibration analysis is reported. These gaps are load-bearing for the claim that the uncertainty estimates 'track' reconstruction errors, and they can be addressed within the manuscript's scope by adding pointwise calibration diagnostics and by making the model-selection procedure more transparent.","major_comments":[{"comment":"The central claim that hierarchical GP uncertainty 'consistently adapted to the amount and quality of the available simulation data' is supported only by visual comparison of averaged RMSE heatmaps and averaged predictive standard deviation heatmaps. No pointwise calibration, empirical coverage, or standardized-residual analysis is reported. Because the method outputs a full predictive covariance, the paper should quantify, for example, the frequency with which the WHAM reference falls inside the 68% or 95% credible intervals as a function of CV position and data-ablation condition, or report z-scores (f_i - μ_i)/σ_i. Without such diagnostics, the observed correlation between global average SD and average RMSE does not establish that the intervals are locally trustworthy, especially near barriers where GP smoothness assumptions may fail.","section":"§3.2, Fig. 3"},{"comment":"The metadynamics support for the headline claim is a single Pearson correlation coefficient ρ=0.95 between average predictive standard deviation and average RMSE, computed over a small number of trajectory-length aggregates (the text lists 1% and then 'remaining trajectory lengths', but the figure appears to show five points). This is not a sufficient quantitative basis for the claim. The correlation has no associated uncertainty or significance test, it averages over all CV positions, and it does not assess whether the predictive intervals have correct coverage. The authors should report the number of points, a pointwise analysis, and preferably a calibration curve or coverage statistic.","section":"§3.4, Fig. 5"},{"comment":"There is a circularity concern between model selection and evaluation. The LOO pseudo-likelihood was selected over the standard marginal likelihood (SI Section D) because it 'more consistently tracked the observed accuracy' on the same RMSE-vs-SD ablation grids used to evaluate the method, and the metadynamics discretization parameters n_h=120 and n_d=40 were chosen partly to keep inferred noise hyperparameters 'of comparable magnitude' and to avoid unstable reconstructions on those same benchmarks. This does not invalidate the method, but it means the reported 'tracking' is partly a result of tuning the objective and preprocessing to the benchmark. To support the general claim, the authors should either report the selection procedure explicitly as exploratory, or validate the method on an independent system with pre-specified settings.","section":"§2.4, SI Sections D and F.2"},{"comment":"The stationary squared-exponential kernel and scalar inferred noise hyperparameters are acknowledged in the Discussion to be a potential source of model mismatch: 'there is no physical reason why a free energy profile should exhibit uniform correlation structure throughout an entire CV domain.' This is more than a conceptual limitation; it directly threatens the UQ claim. If local length scales vary near barriers or basins, the inferred scalar noise and length scale can absorb the mismatch, producing intervals that are locally over- or under-confident even if the average SD correlates with average RMSE. The paper should test this sensitivity, for example by comparing with a non-stationary kernel on at least one benchmark, or by reporting pointwise coverage in data-rich regions where hyperparameters are well identified. The current aggregate evaluation cannot detect such local failures.","section":"§2.3, Eq. (23); Discussion, §4"}],"minor_comments":[{"comment":"Typo: 'mulistate Bennet acceptance ratio' should be 'multistate Bennett acceptance ratio'.","section":"§1, Introduction"},{"comment":"The notation P_w(ξ_i) is used both for a histogram estimate and, later, for probability in the SI noise formulas. Clarifying this would avoid confusion.","section":"§2.1, Eq. (2)"},{"comment":"The table caption says 'LOO MAP hyperparameter estimates' but the text refers to 'generalized hyperposterior MAP'. Please make the terminology consistent.","section":"§3.1, Table 1"},{"comment":"Typo: 'hyperparamter' should be 'hyperparameter'.","section":"SI Section E"},{"comment":"The heatmap color bars are on a logarithmic scale and the white threshold is defined as RMSE=5 kJ/mol / SD=5 kJ/mol. It would be helpful to state explicitly in the caption that the white regions denote values below the threshold, to avoid ambiguity.","section":"§3.2, Fig. 3"},{"comment":"The GitHub and Zenodo links are a welcome addition. Please state the software version used (e.g., Pyro version) and any relevant dependencies, or point to an environment file, to improve reproducibility.","section":"§7, Data Availability"}],"recommendation":"major_revision","confidential_remarks":"The manuscript presents a well-constructed hierarchical GP framework, and the HMC diagnostics and sensitivity analyses are exemplary. The main issue is that the central validation is not yet quantitative enough: the umbrella-sampling evidence is visual, the metadynamics evidence is a single aggregate correlation, and no pointwise coverage or calibration analysis is reported. The authors should be asked to add such diagnostics, to report the number of points underlying the metadynamics correlation, and to address the circularity between model selection and evaluation. If these are addressed, the paper could be suitable for publication. In its current form, the claim that hierarchical GP inference is 'the only approach' with reliably adaptive uncertainty is too strong relative to the evidence presented."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Read this one if you care about convergence judgment in enhanced-sampling workflows. The paper does something useful: it builds a fully Bayesian GP for one-dimensional free-energy profiles, treating observation noise as an inferred hyperparameter, using joint function and derivative observations, and marginalizing kernel hyperparameters via HMC-NUTS. The central empirical claim—that hyperparameter-marginalized GP uncertainty adapts to data quantity while fixed-hyperparameter GP and umbrella integration do not—is new in the GP-free-energy literature, and the authors support it with code, data, and careful MCMC diagnostics. The noise-as-inferred-hyperparameter fix for pathological length-scale modes is the most convincing part; Figure 2 is a clear demonstration.\n\nCredit where due: hyperprior sensitivity, comparison of LOO vs marginal likelihood, MAP vs full sampling, and the metadynamics discretization analysis are all honest pieces of work. The paper is transparent about the stationary-kernel assumption, which matters.\n\nThe soft spot is the main validation. The ablation study reports averaged RMSE and averaged predictive SD heatmaps, but no pointwise coverage, calibration, or standardized-residual analysis. So the strong statement in Section 5—that this was the only approach whose uncertainties 'consistently adapted'—is not actually demonstrated. A method can have average SD tracking average RMSE while being locally over- or under-confident; the acknowledged stationary kernel makes that a real possibility near barriers. The metadynamics support is a single Pearson rho on a small number of trajectory-length aggregates, which is weak. The LOO objective was chosen on the basis of the same benchmark (the authors are upfront about this in SI Section D), and the metadynamics binning parameters were selected using the data—not fatal, but it tightens the inference.\n\nNone of this is a load-bearing flaw. The method is sensible, the math is standard, and the practical value for automated workflows is real. What's missing is a quantitative check of the UQ claim itself: report pointwise coverage of the credible intervals against true errors, or standardized residuals, across the ablation grid. That would make the headline claim defensible.\n\nThis paper deserves a serious referee. I'd send it out, with a request for a calibration analysis and a more tempered Section 5 claim until it exists. I would cite it for the method, not yet for the tracking claim.","headline":"Hierarchical GP free-energy UQ that propagates hyperparameter uncertainty is a real practical step forward, but the signature claim—uncertainty tracks error—rests on visual heatmaps and needs calibration-style evidence before it should be cited as established.","tokens_in":22263,"tokens_out":2323,"would_cite":true,"duration_ms":21125,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"Generalized hierarchical Bayesian inference makes free-energy uncertainty estimates adapt to the amount and quality of simulation data.","keywords":["free energy calculations","uncertainty quantification","Gaussian process regression","hierarchical Bayesian inference","umbrella sampling","metadynamics","hyperparameter uncertainty","peptide-lipid membrane interactions"],"falsifier":"Run repeated independent enhanced-sampling simulations on a system with a known free energy surface (for example a double well with a narrow barrier), apply the same trajectory/window ablation, and compute the empirical coverage of the reported ±1σ credible intervals. If the intervals contain the true profile in substantially fewer than about 68% of reconstructions in the low-data cells, the central claim that the uncertainty tracks real error is falsified.","tokens_in":21439,"feed_emoji":"📊","tokens_out":9699,"duration_ms":75217,"temperature":0.7,"pith_summary":"The paper sets out to show that standard ways of putting error bars on computed free energy profiles — umbrella integration with block averaging, and Gaussian processes with fixed hyperparameters — produce uncertainty estimates that do not respond to how much simulation data was collected. The proposed fix is a generalized hierarchical Gaussian process — a probability distribution over smooth free energy surfaces — in which the kernel hyperparameters and observation noise are treated as unknowns, given priors, and sampled by Hamiltonian Monte Carlo rather than fixed or optimized to a point. On a 100-point data-ablation grid for a peptide–membrane system and on a metadynamics trajectory scan, the hierarchical model's predictive standard deviation tracks the actual reconstruction error, while the baselines either stay nearly constant or collapse abruptly. This matters because a trustworthy convergence signal is what would let automated free energy workflows stop sampling once the answer is known well enough.","feed_headline":"Sample GP hyperparameters, or error bars stay blind to new data","feed_subtitle":"Unlike standard GP and umbrella integration, hierarchical GP credible intervals grow and shrink with data.","key_machinery":"The central object is a generalized hierarchical Gaussian process for the free energy profile: a zero-mean GP with a squared-exponential kernel (a covariance function controlling smoothness and correlation length), fed by both histogram free-energy observations and derivative (mean-force) observations in a joint function–derivative observation model. The kernel length scale, kernel amplitude, and the function and derivative noise variances are given log-normal hyperpriors and sampled from a generalized hyperposterior using Hamiltonian Monte Carlo; the likelihood is a leave-one-out pseudo-likelihood. Hyperparameter samples are propagated into the predictive distribution through the law of tot","core_discovery":"The central claim is that predictive uncertainty in Gaussian-process free energy reconstruction is trustworthy only when the sources of model uncertainty — the GP hyperparameters and the observation noise — are marginalized over rather than fixed in advance or collapsed to a maximum-a-posteriori point. The evidence is a data-ablation study on the R9 peptide–membrane free energy profile and a metadynamics benchmark on phenol–membrane permeation. In both settings the hierarchical GP's predictive standard deviation evolves with the amount and quality of the simulation data and tracks the reconstruction RMSE (Pearson correlation 0.95 in the metadynamics case), whereas a fixed-hyperparameter GP g","pith_inferences":["If the demonstrated convergence signal holds beyond these two systems, the method could serve as a stopping rule in automated adaptive sampling, halting simulations once credible-interval widths fall below a target — a step the paper motivates but does not itself implement.","The stationary squared-exponential kernel is the premise most likely to fail first: real free energy surfaces may have different correlation lengths near barriers than in basins, and a non-stationary kernel comparison on a profile with known narrow barriers would reveal whether inferred noise absorbs the mismatch.","The 'only approach whose uncertainty adapted' claim is made relative to the three baselines tested; a sharper test would pit the hierarchical GP against other adaptive UQ schemes on the same ablation grid.","The leave-one-out-versus-marginal-likelihood finding suggests that for uncertainty-focused tasks, predictive objectives should be preferred over fit-based objectives; this lesson plausibly transfers to other Gaussian-process applications in computational chemistry."],"forward_implications":["A hyperparameter-sampled GP provides a usable convergence signal: its credible intervals narrow as trajectory length and window count increase, letting users distinguish a converged profile from an under-sampled one.","Fixed-hyperparameter GP uncertainty stays essentially constant across the ablation grid, so it cannot indicate whether additional sampling improved the reconstruction.","Umbrella integration with block averaging is underconfident in moderate-data regimes and collapses abruptly once data are nearly complete, making it a misleading convergence diagnostic.","Optimizing hyperparameters rather than sampling them underdisperses the predictive uncertainty by about 18% on average in the low-data setting, and full hyperposterior propagation corrects the overconfidence.","The framework transfers from umbrella sampling to extended Lagrangian metadynamics, where the hierarchical GP's average uncertainty tracks RMSE (correlation 0.95), suggesting it applies to enhanced-sampling methods generally."],"fun_headline_variants":["Marginalize GP hyperparameters or your error bars stay blind","Hierarchical GP error bars track data, fixed ones don't","Sample hyperparameters for GP error bars that adapt to data","Trustworthy free energy errors need sampled hyperparameters"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"The free energy profile is modeled as a zero-mean Gaussian process with a stationary squared-exponential covariance, meaning a single correlation length applies across the whole collective-variable domain; if the true profile has different smoothness around barriers versus basins, the inferred noise will absorb the mismatch and the reported uncertainties may not track true error.","fun_headline_variants_meta":{"raw":{"variants":["Marginalize GP hyperparameters or your error bars stay blind","Hierarchical GP error bars track data, fixed ones don't","Sample hyperparameters for GP error bars that adapt to data","Trustworthy free energy errors need sampled hyperparameters"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000165,"raw_usage":{"total_tokens":1038,"prompt_tokens":650,"completion_tokens":388,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":394,"completion_tokens_details":{"reasoning_tokens":321}},"tokens_in":394,"tokens_out":388,"duration_ms":4343,"temperature":1.0,"reasoning_tokens":321,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-01T05:03:42.572402+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run repeated independent enhanced-sampling simulations on a system with a known free energy surface (for example a double well with a narrow barrier), apply the same trajectory/window ablation, and compute the empirical coverage of the reported ±1σ credible intervals. If the intervals contain the true profile in substantially fewer than about 68% of reconstructions in the low-data cells, the central claim that the uncertainty tracks real error is falsified.","supporting_citations":[],"review_version":1}