Exploring cocoa bean fermentation mechanisms by kinetic modelling

Compared with other fermentation processes in food industry, cocoa bean fermentation is uncontrolled and not standardized. A detailed mechanistic understanding can therefore be relevant for cocoa bean quality control. Starting from an existing mathematical model of cocoa bean fermentation we analyse five additional biochemical mechanisms derived from the literature. These mechanisms, when added to the baseline model either in isolation or in combination, were evaluated in terms of their capacity to describe experimental data. In total, we evaluated 32 model variants on 23 fermentation datasets. We interpret the results from two perspectives: (1) success of the potential mechanism, (2) discrimination of fermentation protocols based on estimated parameters. The former provides insight in the fermentation process itself. The latter opens an avenue towards reverse-engineering empirical conditions from model parameters. We find support for two mechanisms debated in the literature: consumption of fructose by lactic acid bacteria and production of acetic acid by yeast. Furthermore, we provide evidence that model parameters are sensitive to differences in the cultivar, temperature control and usage of steel tanks compared with wooden boxes. Our results show that mathematical modelling can provide an alternative to standard chemical fingerprinting in the interpretation of fermentation data.

If you have been asked to revise the written English in your submission as a condition of publication, you must do so, and you are expected to provide evidence that you have received language editing support. The journal would prefer that you use a professional language editing service and provide a certificate of editing, but a signed letter from a colleague who is a native speaker of English is acceptable. Note the journal has arranged a number of discounts for authors using professional language editing services (https://royalsociety.org/journals/authors/benefits/language-editing/).

===PREPARING YOUR REVISION IN SCHOLARONE===
To revise your manuscript, log into https://mc.manuscriptcentral.com/rsos and enter your Author Centre -this may be accessed by clicking on "Author" in the dark toolbar at the top of the page (just below the journal name). You will find your manuscript listed under "Manuscripts with Decisions". Under "Actions", click on "Create a Revision". Attach your point-by-point response to referees and Editors at Step 1 'View and respond to decision letter'. This document should be uploaded in an editable file type (.doc or .docx are preferred). This is essential.
Please ensure that you include a summary of your paper at Step 2 'Type, Title, & Abstract'. This should be no more than 100 words to explain to a non-scientific audience the key findings of your research. This will be included in a weekly highlights email circulated by the Royal Society press office to national UK, international, and scientific news outlets to promote your work.

At
Step 3 'File upload' you should include the following files: --Your revised manuscript in editable file format (.doc, .docx, or .tex preferred). You should upload two versions: 1) One version identifying all the changes that have been made (for instance, in coloured highlight, in bold text, or tracked changes); 2) A 'clean' version of the new manuscript that incorporates the changes made, but does not highlight them. --If you are requesting a discretionary waiver for the article processing charge, the waiver form must be included at this step.
--If you are providing image files for potential cover images, please upload these at this step, and inform the editorial office you have done so. You must hold the copyright to any image provided.
--A copy of your point-by-point response to referees and Editors. This will expedite the preparation of your proof.

At
Step 6 'Details & comments', you should review and respond to the queries on the electronic submission form. In particular, we would ask that you do the following: --Ensure that your data access statement meets the requirements at https://royalsociety.org/journals/authors/author-guidelines/#data. You should ensure that you cite the dataset in your reference list. If you have deposited data etc in the Dryad repository, please include both the 'For publication' link and 'For review' link at this stage.
--If you are requesting an article processing charge waiver, you must select the relevant waiver option (if requesting a discretionary waiver, the form should have been uploaded at Step 3 'File upload' above).
--If you have uploaded ESM files, please ensure you follow the guidance at https://royalsociety.org/journals/authors/author-guidelines/#supplementary-material to include a suitable title and informative caption. An example of appropriate titling and captioning may be found at https://figshare.com/articles/Table_S2_from_Is_there_a_trade-off_between_peak_performance_and_performance_breadth_across_temperatures_for_aerobic_sc ope_in_teleost_fishes_/3843624.

At
Step 7 'Review & submit', you must view the PDF proof of the manuscript before you will be able to submit the revision. Note: if any parts of the electronic submission form have not been completed, these will be noted by red message boxes.

RSOS-210274.R1 (Revision)
Review form: Reviewer 1 Is the manuscript scientifically sound in its present form? Yes

Are the interpretations and conclusions justified by the results? Yes
Is the language acceptable? Yes

Recommendation?
Accept with minor revision (please list in comments)

Comments to the Author(s)
This revised version of the paper addresses correctly the concerns of the original review. In my opinion, the paper is ready for publication. I have a few additional minor comments: -There is still a problem with considering that "few thousand iterations" is sufficient for "convergence". Note that even with independent MC sampling, one needs to take into consideration the variance of the statistic being approximated. The use of magic numbers like 1,000 or 100,000 samples is widespread but is simply wrong.
-I'm ok with the authors a practical approach to defining successful fit for which the MCMC-NUTS did not fail. With the new wording, this seems reasonable. However, please clarify that this is due to the limitations of the MCMC-NUTS: The authors of Stan would like to make us believe that if their sampler cannot handle a posterior distribution, then the user's posterior is wrong and the user most change their priors, models etc. This is a ridiculous and very convenient argument by the authors of Stan to justify their questionable sampler (which the authors of this paper are all very keen to defend).
-It will be very illustrating to now mention in the paper that, regarding the PCA analysis, inclusion of the uncertainty in posterior, using the MCMC draws, does make a difference in the analysis.

Decision letter (RSOS-210274.R1)
We hope you are keeping well at this difficult and unusual time. We continue to value your support of the journal in these challenging circumstances. If Royal Society Open Science can assist you at all, please don't hesitate to let us know at the email address below.

Dear Mr Moreno-Zambrano
On behalf of the Editors, we are pleased to inform you that your Manuscript RSOS-210274.R1 "Exploring cocoa bean fermentation mechanisms by kinetic modelling" has been accepted for publication in Royal Society Open Science subject to minor revision in accordance with the referees' reports. Please find the referees' comments along with any feedback from the Editors below my signature.
We invite you to respond to the comments and revise your manuscript. Below the referees' and Editors' comments (where applicable) we provide additional requirements. Final acceptance of your manuscript is dependent on these requirements being met. We provide guidance below to help you prepare your revision.
Please submit your revised manuscript and required files (see below) no later than Friday 14 January. Note: the ScholarOne system will 'lock' if submission of the revision is attempted after the deadline. If you do not think you will be able to meet this deadline please contact the editorial office immediately.
Please note article processing charges apply to papers accepted for publication in Royal Society Open Science (https://royalsocietypublishing.org/rsos/charges). Charges will also apply to papers transferred to the journal from other Royal Society Publishing journals, as well as papers submitted as part of our collaboration with the Royal Society of Chemistry (https://royalsocietypublishing.org/rsos/chemistry). Fee waivers are available but must be requested when you submit your revision (https://royalsocietypublishing.org/rsos/waivers). -There is still a problem with considering that "few thousand iterations" is sufficient for "convergence". Note that even with independent MC sampling, one needs to take into consideration the variance of the statistic being approximated. The use of magic numbers like 1,000 or 100,000 samples is widespread but is simply wrong.
-I'm ok with the authors a practical approach to defining successful fit for which the MCMC-NUTS did not fail. With the new wording, this seems reasonable. However, please clarify that this is due to the limitations of the MCMC-NUTS: The authors of Stan would like to make us believe that if their sampler cannot handle a posterior distribution, then the user's posterior is wrong and the user most change their priors, models etc. This is a ridiculous and very convenient argument by the authors of Stan to justify their questionable sampler (which the authors of this paper are all very keen to defend).
-It will be very illustrating to now mention in the paper that, regarding the PCA analysis, inclusion of the uncertainty in posterior, using the MCMC draws, does make a difference in the analysis.

===PREPARING YOUR MANUSCRIPT===
Your revised paper should include the changes requested by the referees and Editors of your manuscript.
You should provide two versions of this manuscript and both versions must be provided in an editable format: one version should clearly identify all the changes that have been made (for instance, in coloured highlight, in bold text, or tracked changes); a 'clean' version of the new manuscript that incorporates the changes made, but does not highlight them. This version will be used for typesetting.
Please ensure that any equations included in the paper are editable text and not embedded images.
Please ensure that you include an acknowledgements' section before your reference list/bibliography. This should acknowledge anyone who assisted with your work, but does not qualify as an author per the guidelines at https://royalsociety.org/journals/ethicspolicies/openness/.
While not essential, it will speed up the preparation of your manuscript proof if you format your references/bibliography in Vancouver style (please see https://royalsociety.org/journals/authors/author-guidelines/#formatting). You should include DOIs for as many of the references as possible.
If you have been asked to revise the written English in your submission as a condition of publication, you must do so, and you are expected to provide evidence that you have received language editing support. The journal would prefer that you use a professional language editing service and provide a certificate of editing, but a signed letter from a colleague who is a proficient user of English is acceptable. Note the journal has arranged a number of discounts for authors using professional language editing services (<a href="https://royalsociety.org/journals/authors/benefits/languageediting/">https://royalsociety.org/journals/authors/benefits/language-editing/</a>).

===PREPARING YOUR REVISION IN SCHOLARONE===
To revise your manuscript, log into https://mc.manuscriptcentral.com/rsos and enter your Author Centre -this may be accessed by clicking on "Author" in the dark toolbar at the top of the page (just below the journal name). You will find your manuscript listed under "Manuscripts with Decisions". Under "Actions", click on "Create a Revision".
Attach your point-by-point response to referees and Editors at the 'View and respond to decision letter' step. This document should be uploaded in an editable file type (.doc or .docx are preferred). This is essential, and your manuscript will be returned to you if you do not provide it.
Please ensure that you include a summary of your paper at the 'Type, Title, & Abstract' step. This should be no more than 100 words to explain to a non-scientific audience the key findings of your research. This will be included in a weekly highlights email circulated by the Royal Society press office to national UK, international, and scientific news outlets to promote your work. An effective summary can substantially increase the readership of your paper.
At the 'File upload' step you should include the following files: --Your revised manuscript in editable file format (.doc, .docx, or .tex preferred). You should upload two versions: 1) One version identifying all the changes that have been made (for instance, in coloured highlight, in bold text, or tracked changes); 2) A 'clean' version of the new manuscript that incorporates the changes made, but does not highlight them. --If you are requesting a discretionary waiver for the article processing charge, the waiver form must be included at this step.
--If you are providing image files for potential cover images, please upload these at this step, and inform the editorial office you have done so. You must hold the copyright to any image provided.
--A copy of your point-by-point response to referees and Editors. This will expedite the preparation of your proof.
At the 'Details & comments' step, you should review and respond to the queries on the electronic submission form. In particular, we would ask that you do the following: --Ensure that your data access statement meets the requirements at https://royalsociety.org/journals/authors/author-guidelines/#data. You should ensure that you cite the dataset in your reference list. If you have deposited data etc in the Dryad repository, please only include the 'For publication' link at this stage. You should remove the 'For review' link.
--If you are requesting an article processing charge waiver, you must select the relevant waiver option (if requesting a discretionary waiver, the form should have been uploaded, see 'File upload' above).
--If you have uploaded any electronic supplementary (ESM) files, please ensure you follow the guidance at https://royalsociety.org/journals/authors/author-guidelines/#supplementarymaterial to include a suitable title and informative caption. An example of appropriate titling and captioning may be found at https://figshare.com/articles/Table_S2_from_Is_there_a_trade-off_between_peak_performance_and_performance_breadth_across_temperatures_for_aerobic_sc ope_in_teleost_fishes_/3843624. At the 'Review & submit' step, you must view the PDF proof of the manuscript before you will be able to submit the revision. Note: if any parts of the electronic submission form have not been completed, these will be noted by red message boxes -you will need to resolve these errors before you can submit the revision.

Decision letter (RSOS-210274.R2)
We hope you are keeping well at this difficult and unusual time. We continue to value your support of the journal in these challenging circumstances. If Royal Society Open Science can assist you at all, please don't hesitate to let us know at the email address below.
Dear Mr Moreno-Zambrano, I am pleased to inform you that your manuscript entitled "Exploring cocoa bean fermentation mechanisms by kinetic modelling" is now accepted for publication in Royal Society Open Science.
Please remember to make any data sets or code libraries 'live' prior to publication, and update any links as needed when you receive a proof to check -for instance, from a private 'for review' URL to a publicly accessible 'for publication' URL. It is good practice to also add data sets, code and other digital materials to your reference list.
Our payments team will be in touch shortly if you are required to pay a fee for the publication of the paper (if you have any queries regarding fees, please see https://royalsocietypublishing.org/rsos/charges or contact authorfees@royalsociety.org).
The proof of your paper will be available for review using the Royal Society online proofing system and you will receive details of how to access this in the near future from our production office (openscience_proofs@royalsociety.org). We aim to maintain rapid times to publication after acceptance of your manuscript and we would ask you to please contact both the production office and editorial office if you are likely to be away from e-mail contact to minimise delays to publication. If you are going to be away, please nominate a co-author (if available) to manage the proofing process, and ensure they are copied into your email to the journal.
Please see the Royal Society Publishing guidance on how you may share your accepted author manuscript at https://royalsociety.org/journals/ethics-policies/media-embargo/. After publication, some additional ways to effectively promote your article can also be found here https://royalsociety.org/blog/2020/07/promoting-your-latest-paper-and-tracking-yourresults/.
On behalf of the Editors of Royal Society Open Science, thank you for your support of the journal and we look forward to your continued contributions to Royal Society Open Science. The authors present a comparative study of a base line ODE model for cocoa fermentation with several other models varying in complexity and qualitative differences. The adding complexities in the mechanistic model are nicely summarized in Figure 1, with 5 added mechanisms into the fermentation process signifying adding terms to the rhs of the system of ODE's. The 5 added variants are consider to create alternative models, leading to 31 possible alternatives along with the base line model.
1. L. 78-83: Note that proper data modeling and careful consideration of parameter units are common strategies for obtaining 'comparable' estimates from different data sets. Data pre-procesing may well be sound, but the authors need to explain in more detail why their models do not report inherent 'comparable' parameters and why this preprocessing is preferred and correct. Can the models be re-parametrized to have comparable parameters? ss/1009212519 for example. All other approaches are approximations to the latter (e.g. the BIC, DIC etc.). The problem is estimating such normalizations constants using MCMC, which remains a difficult problem, still, see https://doi.org/10.1198/106186008X320258 . The authors should try to use this formal approach, before resorting to other more ad hoc or heuristic methods. On the contrary, the posterior may be very simple (nice unimodal bell shape) to sample from, the MCMC be excellent but the model itself be useless.
Conceptually, the authors need to separate MCMC convergence and model fit. The former needs to be assured first, if no convergence is attained, a larger sample size may be needed or a different sampler altogether. There are other multipurpose samplers that may be used, like DRAM or the t-walk MCMC samplers. Once MCMC convergence has been attained then model adequacy may be analyzed.
I am confident that the authors results will hold, in general, since a detailed process of model assessment was done, but a more methodologically and conceptually sound analysis is in order.

Exploring cocoa bean fermentation mechanisms by kinetic modelling
Mauricio Moreno-Zambrano 1 , Matthias S. Ullrich 1 , and Marc-Thorsten Hütt 1 The authors present a comparative study of a base line ODE model for cocoa fermentation with several other models varying in complexity and qualitative differences. The adding complexities in the mechanistic model are nicely summarized in Figure 1, with 5 added mechanisms into the fermentation process signifying adding terms to the rhs of the system of ODE's. The 5 added variants are consider to create alternative models, leading to 31 possible alternatives along with the base line model.
1. L. 78-83: Note that proper data modeling and careful consideration of parameter units are common strategies for obtaining 'comparable' estimates from different data sets. Data pre-procesing may well be sound, but the authors need to explain in more detail why their models do not report inherent 'comparable' parameters and why this preprocessing is preferred and correct.
We agree with the reviewer that our data-processing steps need more clarification. Answering part by part their concerns, we would like to mention that we indeed reported comparable parameters with those available for single-strained microbial growth studies in literature. Specifically, we made use of a re-scaled version of the obtained parameters of model iteration (MI) involving mechanisms 2 and 3 (MI(2,3)) in Section 4.2 and Supplementary Table S6 (now Supplementary Table S8, after reviewer's comments) when discussing parameter agreement between those obtained from our Bayesian framework and their corresponding values for singlestrained microbial experiments. Thus, in order to give the reader a better explanation, we have changed paragraph mentioned by the reviewer as follows: In all cases, population growth of Y, LAB and AAB were transformed from log base 10 of colony forming units (log 10 (CFU)) to milligrams of microbial group (MG) per gram of pulp (mg(MG) g(pulp) −1 ) as these are the units in which most of kinetic single-strained microbial growth studies report their dynamics as well as their dependent constants, i.e., specific maximum growth and mortality rates, and yield coefficients. Moreover, with the purpose of facilitating the estimation of models' parameters by avoiding numerical issues during their calibration, all time series were scaled by dividing each observation by its own maximum value. Hence, obtained parameter estimates were re-scaled to their original units by using simple transformations for their further comparison with previously reported values for single-strained microbial studies (see Supplementary Material, Section 1) [1].
Furthermore, we have added in the Supplementary material an extra section (now as Section 1 of this material) where we illustrate a simple example of how the re-scaling of the obtained parameters is performed to their original units as follows:

Parameter estimates re-scaling to their original units
Consider the kinetics of a given x variable in any of our proposed models, its concentration [x] at any time and its maximum concentration [x] max within the whole interval in which the fermentation took place.
We can then let x ′ represent a transformation of x as then, the time derivative of x ′ will be given by where [x] max is the maximum concentration of any of the state variables depending on which ODE of a proposed model is expressed in the form of Eq. (2).
Generalizing, one can express all ODEs of a model in the form of Eq. (2). As an example, let's express Eq. (1.MM) for the time derivative of glucose (Glc) as a scaled version in accordance to Eq. (2) where µ 1 , µ 3 , Ks 1 , Ks 3 , Y 1 and Y 2 are the solutions of this ODE parameters upon the scaled time series of Glc (according to Eq. (1)) for specific maximum growth rates (µ i ), substrate saturation constants (Ks i ) and yield coefficients (Y i ) respectively.
In other words, all aforementioned terms are scaled posterior distributions resulting from our Bayesian framework for which re-scaled versions to their original units can be accomplished by properly working out the scaled ODE.
Hence, Eq. 3 can be simplified to where the terms µ 1 , µ 3 , Similarly, transformation factors as those previously described can be then obtained for the rest of parameters in any given MI presented in this study.
Can the models be re-parametrized to have comparable parameters? Scaling of the data obtained from each of the fermentation trials represents itself a re-parametrization of the models in the way we have set up the Bayesian framework. In this sense, due to the differences between orders of magnitude of time series reported for microbial growth and metabolites (e.g., the maximum concentrations of acetic acid bacteria and ethanol by Camu et al. [2] are 0.0019 mg g(pulp) −1 and 22.4920 mg g(pulp) −1 respectively [1]), data scaling offers the advantage of setting a more general scheme of choosing priors within a Bayesian framework. In this way we do not require other model re-parametrizations that could end up with an increasing number of parameters (due to individual priors for each variable). For example, a usual practice of a simpler re-parametrization for parameters describing the dynamics of state variables, which highly differ between each other, can be done by multiplying microbial-related parameters (e.g., yield coefficients and substrate saturation constants) by arbitrary values that could lesser the effect of the magnitude in which metabolites differ from the measurements of microbial growth. However, this approach becomes difficult to implement when fitting several datasets, where those arbitrary values need to be set one by one depending on how much do the microbial time series differ with respect to their corresponding metabolite kinetics across different trials leading also to more complicated prior values choices. Finally, even if such an approach would succeed, will also require a transformation to the parameter real units for the sake of comparison. The latter, as we have above described, is also possible with our re-parametrization choice.
2. Diffuse priors may be become a big problem in model comparison; extensive literature has been written by the Bayesian community on this respect. Please explain prior selection in more detail. We agree with the reviewer that using diffuse priors may become a problem in model comparison. However, the limited knowledge regarding the kinetics of the cocoa bean fermentation process impedes the use of more informative priors. Besides, we consider that prior distributions in a regularized parameter space as the depicted in Section 2.4.2 (thanks to the scaling of the original data) have a weakly informative nature rather than a diffuse one. Considering this, and following the reviewer's concern, we have rephrased Section 2.4.2 as follows.

Choice of priors
The regularization procedure of the data described in Section 2.1 permits to reduce the parameter search space in a convenient way for choosing the priors in ranges between 0 to ≈1, which brings three main advantages: (1) independent prior distributions for each parameter can take the same form, (2) by introducing scale information of the original units in which the parameters of the models are originally measured, we can formulate weakly informative priors capable of covering all possible values in the scaled space [1,3], and (3) by avoiding diffuse priors, further model comparisons will be less likely of being affected by common problems, such as over-fitting [4] and ill-defined posteriors [5].
Hence, posterior distributions of θ and σ were computed using a normal distribution with mean 0.5 and standard deviation of 0.3 for each element of the parameter vector θ and a Cauchy distribution with location 0 and scale of 1 for σ. With the purpose of avoiding estimates with negative values, both priors were truncated to the positive set of real numbers, For a detailed description of prior distributions re-scaled to the parameters' original units see Supplementary Table S2. 3. Why RK45 is used. What time step size was used? This is a very old method, far more advance methods are available with implicit and automatic time steps. We agree with the reviewer that newer and more advance methods are available (e.g., EK1 and EK0 [6]). However, implementation of new ODE solvers in Stan is not a straightforward task because of the need of their translation to Stan language [7]. Besides, we would like to emphasize that Stan makes use of the No-U-Turn Sampler (NUTS) to sample from a target distribution. As NUTS is based on the Hamiltonian Monte Carlo (HMC) algorithm [8]. Among other requirements, this implies a need of computing the gradient of the logposterior by automatic differentiation before starting the sampler [9]. As a result, in the case of implementing an ODE solver with automatic time steps, it will make necessary to update the Jacobian matrix of the ODE model and hence, the log-posterior, at any time there is a change in step size leading to an unavoidable slowdown of computation time and also to a lengthy verification process for testing the adequacy of such a newer ODE solver (e.g., running tests and solving benchmark problems). Instead, we preferred to make use of the built-in Stan solver rk45 that has been extensively used in a wide range of problems and applications (see Rosenbaum et al. [10], Bandiera et al. [11], Grinsztajn et al. [12] to name a few). We do also agree that it is necessary to specify the time step size used, reason why we have added in line 194 of our manuscript updated version, the following: . . . with relative and absolute tolerance values of 1 × 10 −6 for both, and a maximum number of steps of 1 × 10 4 . 4. L. 189: Why the burn in or warm up was left at 1,000 iterations? Why such a small sample size? (4 parallel chains of 3,000, why use parallel chains?). We understand the concern of the reviewer in our apparently low number of iterations used for both, warm-up and sampling of the posterior distributions. Contrary to other sampling algorithms (e.g., Metropolis-Hastings and Gibbs), we here use a sampling, which explores the geometry of the posterior distribution in a more efficient way than purely random walk based approaches. This is due to its formulation on basis of the HMC algorithm [8,9]. As a result, convergence using NUTS can be often reached after few thousand iterations, even for complex models [13]. About the use of 4 parallel chains, this decision is based on two main arguments: (1) running in parallel multiple chains with different starting values allows to check for convergence of independent chains once stationarity is reached (as it is usually represented by traceplots where the called "caterpillar" shape is accomplished) [14]. And, (2) four is the minimum number of multiple chains suggested for using MCMC-NUTS. This responds to the need of counting with multiple chains for diagnostic statistics (i.e., R, tail effective sample size and bulk effective sample size) that perform better with multiple chains due to their nature of taking into account within and between chain variability [14,15].
Note that this would be quite short for an independent sampler, let alone with the correlation inherent in the MCMC. What is then the effective sample size? Another consequence of using NUTS is a considerable reduction in the correlation between drawn samples from the posterior distribution. Opposed to MCMC samplers based on random walks, due to its Hamiltonian dynamics, NUTS offers a more efficient exploration of the target distribution even in the presence of local correlations where the former samplers tend to get stuck [14,16]. In consequence, adequate effective sample sizes (ESS) are obtained with less iterations. Thus, in response to the reviewer concern, we have added an extra table in the Supplementary Material as follows:  -  - Where the averaged ratio of ESS with the total number of draws (ESS/N draws ) from the posteriors of each MI is presented. As a rule of thumb, it has been suggested that an advisable ESS should be in the thousands [15]. Having in mind that our approach used as basis 8000 iterations in total per each MI, in all successful MIs, ESS is acceptable. Consequently, we added at the end of Section 2.4.3 the following sentence: Assessing autocorrelation of the sampled parameters was performed by means of averaging computed effective sample sizes over number of posterior draws (ESS/N draws ) and checking whether these were above 12.5%, meaning that as minimum 1000 ESS were obtained [15].
Finally, we added the following sentence in line 267 of the updated version of our manuscript, section 3.1: Furthermore, in terms of autocorrelation, all successful fits showed averaged values of ESS/N draws over 12.5% (see Supplementary Table S6), indicating an acceptable ESS above 1000.
5. Note that the corner stone of the paper is the model fit-model assessment. Section 2.5 is only very briefly explained and unclear. The authors should explain in far more detail their approach. For example, what is a "succesful fit" (l. 199)?
We agree with the reviewer that we needed to be more detailed in defining what is a "succesful fit". With this in mind, we have added the following paragraph in Section 2.5: Under the umbrella of the proposed Bayesian framework, it is important to define what will be called from now on a successful fit. For this aspect, we consider as such, any MI fit that converged according to the criteria described in Section 2.4.3 and does not involve any re-parametrization or use of distinct priors in cases were divergences of MCMC-NUTS or complete lack of sampling could arise as result of complicated geometries imposed to the posterior distributions by inclusion of the assessed mechanisms over particular datasets.
Note that there is a formal Bayesian approach to model adequacy-model selection: The posterior probability of each model is proportional to its normalization constant, see https://doi.org/10.1214/ss/1009212519 for example. All other approaches are approximations to the latter (e.g. the BIC, DIC etc.). The problem is estimating such normalizations constants using MCMC, which remains a difficult problem, still, see https://doi.org/10.1198/106186008X320258 . The authors should try to use this formal approach, before resorting to other more ad hoc or heuristic methods. We have followed the reviewer's suggestion of implementing Bayesian Model Averaging (BMA) [17]. For this purpose, we used the weights for each MI that resulted from performing pseudo-BMA [18] over sets of MIs that successfully fitted a given dataset and averaged them over the total number of datasets that were fitted at least once. Hence, we have modified Section 2.5 as follows: On the one hand, Bayesian Model Averaging (BMA [17]) weights were computed using pseudo-BMA [18]. These weights, understood as representations of the relative probability of each MI [19], were then averaged over the total number of datasets that were fitted at least once by any MI, and used as a measure of the adequacy of each MI across all data. For computing the mean value of BMA weights (BMA w ), non-successful fits were assigned values of zero. Furthermore, an observed success rate (OSR) and expected success rate (ESR) were determined on the basis of times where the model was satisfactorily fit to a given dataset. OSR is then defined by the ratio of the number of successful fits and the total of datasets fitted by at least one MI. ESR, used to properly compare the success rates of models with only a single additional mechanism with those models containing a combination of mechanisms, is defined as the product of the OSRs of the elementary MIs. For model variant MI (1,2,5), for example, the expected success rate is then the product of the observed success rates of the elementary models MI(1), MI(2) and MI (5).
As a consequence, we have added in section 3.2 the following paragraphs: Both approaches for assessing model success across all datasets, a formal one as BMA and our proposed measures OSR and ESR, showed similar conclusions over competing MIs.  2,3)). Similarly, among these, MI(2,3) had a maximum OSR of 0.74 as listed in Table 4.
Furthermore, we updated Table 4 in the main manuscript accordingly: Table 4. Summary of successful fits across 31 models iterations (MIs) and baseline. Light green-coloured cells indicate successful fits. Light-red coloured cells indicate non-successful fits. Columns "MI( )", "#", "BMA w ", "OSR" and "ESR" refer to combination of mechanisms deployed in model iteration, number of parameters, averaged Bayesian Model Averaging weights, observed success rate and expected success rate, respectively. The latter is also evident by higher BMA weights by such MIs (see Supplementary Table S7).
And added to the Supplementary Material, Table S7 as follows: Table S7. Summary of Bayesian Model Avering (BMA) among successful fits across 31 model iterations (MI). BMA's weights in the referred MI is reported. Column "MI( )" indicates combination of mechanisms deployed in MI.
We have followed the reviewer recommendation and added a more detailed description of the PCA conduction. Hence, in section 2.6 we have added the following sentence: In an approach of linking our mathematical exploration of mechanisms with a real-life application (where Principal Component Analysis (PCA) is a common approach), we studied, whether projecting the obtained vectors of kinetic parameters on the PCA space allows for a separation (and hence classification) of experiments according to fermentation features (e.g., cocoa beans' country of origin and used cultivar), which are of interest in both, academic research and chocolate industry.
Regarding the reviewers' concern of whether we did use all data sets and all models for performing PCA, we indeed did it, and we have pointed out this by adding the following sentence in Section 2.6: Consequently, six main features of fermentation trials were taken into account as groups and analysed via PCA over parameter estimates for all successful MIs.
Concerning the possibility of using the uncertainty in the posterior for conducting PCAs, we have followed the recommendation of the reviewer. For accomplishing it, we have used all posterior draws within the 95% credible interval (CI), as now stated in Section 2.6 in the following sentence: Only the posterior samples contained within their 95% credible interval (CI) from each chain in the MCMC-NUTS runs were taken into account to perform PCA.
As one could expect, by using the uncertainty of the posteriors for PCAs, led to a change in the reported explained variances that we have updated in the following manner: It is important to mention that conducting PCAs in this way, has not drastically changed our original results and interpretations. We have added minor corrections in light of the updated PCAs as follows:      Figure S28. PCA score (a) and loading plot (b) from all parameters of model iteration MI (2). For visualization purposes, scores of only 10% of posterior draws are shown. Heap (hp), platform (pt), stainless-steel tank (st) and wooden box (wb) fermentation methods are shown. Parameters located on the left and right with respect to 0 in the horizontal axis in the loading plot determine differentiation between stainless-steel tank and rest of fermentation methods, respectively. Finally, we have updated also Figure S27 regarding Mahalanobis distances, which has not changed our original findings either. 7. The authors seem to trust nearly blindly the statistics R etc. The authors are wrong here in believing that a bad convergence in the MCMC means a bad model. The model can be very good but results in a difficult posterior to sample from and the MCMC method itself is what is inadequate. The U-turn sampler is not bullet proof, nor any of the other methods programmed in Stan. Difficult posteriors (multimodal, too correlated) are difficult to sample for any method and the methods in Stan are no exceptions (btw contrary to what the authors of Stan seem to claim). This misconception is also apparent in the Discussion, l. 354. On the contrary, the posterior may be very simple (nice unimodal bell shape) to sample from, the MCMC be excellent but the model itself be useless.
Conceptually, the authors need to separate MCMC convergence and model fit. The former needs to be assured first, if no convergence is attained, a larger sample size may be needed or a different sampler altogether. There are other multipurpose samplers that may be used, like DRAM or the t-walk MCMC samplers. Once MCMC convergence has been attained then model adequacy may be analyzed.
I am confident that the authors results will hold, in general, since a detailed process of model assessment was done, but a more methodologically and conceptually sound analysis is in order.
We agree with the reviewer in that we did not make a clear definition between bad convergence and model fitting. After our corrections, specially the one given as response for comment 5 of the reviewer, we have made clear a definition of what it is considered a successful fit under our Bayesian framework. Taking such a definition as a starting point, we now feel confident that the structure of our manuscript has a well defined final separation between MCMC convergence and model fit by changing "convergence" wording for more appropriate phrases in line with our definition of successful fits in the Discussion section of our manuscript. In more detail, these changes in our updated version are: Line 378: Section 4.1 has been renamed as "Model plausibility". Lines 380, 425, 435: Change of the word "convergence" for "success", "sample from complicated geometries of the posterior target" and "successful fits" respectively. Lines 381, 385, 393, 396, 413, 449: Change of the word "non-convergence" by "failure in fitting" in line 381, and "fitting failure" in the rest of the cases.
Finally, we would like to thank the anonymous reviewer for their advise that has led us to an improved new version of our original manuscript.
In conclusion, we are certain that our chosen number of iterations for fitting each model is in line to what specialized literature suggest as good practices for assessing convergence.
I'm ok with the authors a practical approach to defining successful fit for which the MCMC-NUTS did not fail. With the new wording, this seems reasonable. However, please clarify that this is due to the limitations of the MCMC-NUTS: The authors of Stan would like to make us believe that if their sampler cannot handle a posterior distribution, then the user's posterior is wrong and the user most change their priors, models etc. This is a ridiculous and very convenient argument by the authors of Stan to justify their questionable sampler (which the authors of this paper are all very keen to defend).
We agree with the reviewer that limitations of MCMC-NUTS are also among possible causes for our models to fail. To address this concern, we have added among the possible causes for their failure in line 384 of the Discussion, section 4.1 of the corrected manuscript the following: . . ., and even due to known limitations of MCMC-NUTS, such as poor chain mixing in presence of multimodal posteriors [8] and its higher sensitivity to model parameterization with respect to other samplers [9].
It will be very illustrating to now mention in the paper that, regarding the PCA analysis, inclusion of the uncertainty in posterior, using the MCMC draws, does make a difference in the analysis.
Following the reviewer's recommendation, we have added an extra paragraph in line 506 of the Discussion, section 4.3 of the corrected manuscript as follows: Lastly, it is worth to remark that including all parameters' posterior draws contained in their 95% CI for conducting PCAs plays an important role in finding more meaningful differences between the features. This conclusion comes from a former approach were we performed PCAs over the medians of each chain of successful MIs instead. In such an exercise, besides obvious quantitative differences between its PCs and their counterparts reported here (see Supplementary Table S9), the same classification of the features was observed, with the exception of a noticeable separation between all countries (see Supplementary Figure S31) rather than the overlap between Ghana and Ecuador shown in Supplementary Figure S29, which could have been ignored if it were not by considering the posterior distributions' uncertainties.
Consequently, we have added Supplementary  Finally, we would like to thank again the anonymous reviewer for their advise which had led us to improve this final version of our manuscript.