This is a very timely and impactful paper describing an insightful approach to model-free fitting of relaxation data. I predict that the proposed approach will be hugely influential and likely replace the seminal model-selection protocol published by the Palmer group in 1995, which has some 1000+ citations.

Here Crawley and Palmer implement bootstrap aggregation as a novel way to perform model-free fitting, taking into consideration the specific demands imposed by NMR relaxation data sets, where individual data points (i.e. J(w) sampled at specific w values) have very different impact on the fitted model.

The paper is comprehensive and easy to follow. Well-chosen examples are presented highlighting the advantages of the bootstrap aggregation approach. Nonetheless, I have a few suggestions for slight improvements in presentation and two questions.

1. Could you please comment on the number of field strengths required to successfully implement bootstrap aggregation. Assuming I get the math right, using 3 field strengths, the total number of possible samples is 7³= 343, which seems to be at the low end of the number of samples normally employed, but perhaps this would suffice?
2. How do you determine the correlation time for overall rotational diffusion, tau_m (in the case of a globular protein)? The suggested protocol fits tau_m individually for each residue. Do you forego the concept of fitting a global tau_m as part of the MF fits? And subsequently fit a rotational diffusion model (isotropic or anisotropic) to these individual values (while taking into account the orientation of the HN bond vectors in the molecular frame in the case of anistropic models)?

Minor points:

line 19, Suggestion: spell out Akaike and Bayesian Information Criterion when introducing AIC and BIC.

l. 145-148, the mean J(0) is presumably only used in the conventional MF protocol with MC error analysis(?) This should be stated here to avoid confusion.

The bootstrap aggregation protocol is well described on p. 6, but I still feel that it might be beneficial to include flow-chart type figure outlining the construction of the bootstrap sample datasets.

The tables are not easily interpreted without referring back to the text. Please add footnotes to define p_ij and Y_ij in words (Table 1). Please add text to indicate that “Smooth” refers to the percentage of selected models in Tables 2-4.

Typos: ”paramaters” (l. 86); ”interogating conformional” (l. 250)

Crawley and Palmer propose the application of bootstrap aggregation methodology to derive uncertainty estimates for the analysis of magnetic field strength-dependent NMR relaxation data in the context of the model-free spectral density function selection protocol previously introduced by Palmer and colleagues. A weakness of the least-squares algorithm described in the widely used original protocol is its susceptibility to model-selection error. Sets of relaxation values that are statistically consistent with each other in terms of estimated experimental uncertainty may give rise to the selection of different model-free spectral density functions. In the present manuscript, the statistical method of bootstrap aggregation is used to enable a joint refinement in the parameter uncertainty estimation and the dynamical model selection process.

As applied to ¹⁵N relaxation data for the bZip domain of GCN4, bootstrap aggregation is reported to reduce residue-to-residue variations in optimal model-free parameters, particularly in the partially disordered basic region. Regarding the more general practical utility of the proposed sampling protocol, while experimental relaxation data collected at four magnetic field strengths yields 6859 suitably filtered combinations of bootstrap samples, as noted by the first reviewer, the robustness of the statistical analysis may appreciably decline when this value drops to 343 for three magnetic field strengths, and presumably will decrease significantly further when it drops to only 27 for data from two magnetic field strengths.

A key step in the proposed joint refinement process calculates each of the averaged dynamical parameters by summing over the estimates obtained from each of the five spectral density representations being utilized, as weighted by how often each of these five models have been selected (Eq. 10). A potential concern over this approach arises from the fact that while the same set of symbols (τ_m, S_f²,τ_f,S_s²,τ_s) are utilized in each of the five dynamical models used, the functional significance of each symbol is defined within the context of the specific equation being used. Each of these five model equations that are used to represent the spectral density function is capable of accurately fitting only a small subset of the physically plausible spectral density curves. Systematic bias can potentially arise not only with respect to a given dynamics parameter being utilized in distinct model representations but also as a result of the inadequacy with which each of the five model spectral density equations are capable of representing the physical dynamics of the system. While such biasing effects are surely diminished for Model 4 and 5 which incorporate four and five adjustable parameters, respectively, more promising might be the utilization of alternative model equations for the spectral density function that can more robustly represent the range of motion occurring in protein molecules which utilize a smaller set of adjustable parameters for optimization against experimental relaxation data. 

Excellent manuscript that brings new approaches to the problem of model selection in the analysis of model-free dynamics data. It has been recognized for a long time the problem of selecting the model that most robustly fits the analysis, and then the issue of selecting a single model that is marginally "better". This manuscript uses bootstrap aggregation, from machine learning approaches to improve the possibilities of what is (are) the best model(s) and how that adds to our knowledge of the influsence of dynamics on structure/function - the purpose of doing the work. I particularly liked the Arg examples given - especially Arg11 and also Arg26, the idea that the dynamics reflects its position at the juncture of the coiled-coil and basic regions. More conventional analyses may miss this insight.

Like the other reviewers, discussion needs to address the issues of using two or three fields. Data at two fields are the most commonly acquired and so would this method be inapproriate or do the authors have alternate approaches/ideas. 

The paper is well-written and accessible to most in the field. A good balance of theory, method and application. I agree that the tables need better notes - for example the description of the colour scheme in Figure 1 has to be repeated in Fig 2 and 3; better decriptive footnotes in Table 2,3,4 (and actually I think this could be put into a single table with residue first column and clear breaks); possibly the same for Tables 5,6,7.

Very few errrors found. Including the detected typos line 199 "highlighted"