Medusa: Software to build and analyze ensembles of genome-scale metabolic network reconstructions.

Uncertainty in the structure and parameters of networks is ubiquitous across computational biology. In constraint-based reconstruction and analysis of metabolic networks, this uncertainty is present both during the reconstruction of networks and in simulations performed with them. Here, we present Medusa, a Python package for the generation and analysis of ensembles of genome-scale metabolic network reconstructions. Medusa builds on the COBRApy package for constraint-based reconstruction and analysis by compressing a set of models into a compact ensemble object, providing functions for the generation of ensembles using experimental data, and extending constraint-based analyses to ensemble scale. We demonstrate how Medusa can be used to generate ensembles and perform ensemble simulations, and how machine learning can be used in conjunction with Medusa to guide the curation of genome-scale metabolic network reconstructions. Medusa is available under the permissive MIT license from the Python Packaging Index (https://pypi.org) and from github (https://github.com/opencobra/Medusa), and comprehensive documentation is available at https://medusa.readthedocs.io/en/latest.

This is a PLOS Computational Biology Software paper.

Introduction

Hypothesis-driven computational models of biological systems are being increasingly applied to guide experimentation [1]. In hypothesis-driven modeling, in contrast to data-driven modeling [2], hypothesized biological parts, functions, and interactions are mathematically formalized to allow in silico experimentation. These models take many forms, ranging in complexity from a single linear equation relating two quantities to systems of nonlinear differential equations describing dynamic systems.

Across all hypothesis-driven modeling frameworks, the choice of model scope and parameter values may strongly influence simulation results. For some types of hypothesis-driven models in biology, approaches from other fields have been applied to quantify the influence of parameter values on simulation outcomes, such as sensitivity analysis of dynamical models [3]. For network-based models of biological systems such as metabolic or signaling networks, the presence or absence of a network component may be uncertain due to lack of characterization or uncertainty in data itself. Traditional sensitivity analysis methods have recently been reformulated for these systems to analyze sensitivity to topological variation, but these methods have not seen wide adoption [4]. While uncertainty in network structure poses analytical difficulties, it also presents an actionable framework to accelerate biological discovery. Alternative network structures can guide experimental design, allowing comparison of simulation results for alternative networks to experimental data to identify the network structure most consistent with biological behavior (i.e. model selection) [5]. This uncertainty can also be used to prioritize experiments that will maximally improve confidence in the simulations performed with a model (i.e. uncertainty reduction) [6].

In studies of metabolism, genome-scale metabolic network reconstructions (GENREs) have emerged as a useful formalism for hypothesis-driven modeling [7]. In conjunction with biological objective functions, such as maximization of growth rate, GENREs can be used to construct genome-scale metabolic models (GEMs). In addition to topological uncertainty (e.g., presence/absence of reactions in a network), simulations with GEMs generally yield many alternative solutions. Even the simplest simulations that can be performed with GEMs exhibit this behavior. This is the case for flux balance analysis (FBA), in which a pseudo-steady state is assumed, and flux values are found for all reactions in a GEM such that an objective function is optimized [8]. While a single global maximum value for the objective is guaranteed to be found, flux through every other component of the network is only constrained within a solution space, not to a single value. As a result, even though performing FBA yields a single value for the flux through reactions in a network, there are an infinite number of feasible flux values within the range determined by the solution space for some reactions. Techniques such as flux variability analysis and flux sampling have been developed to explore the space of alternative solutions in this scenario [9,10].

A myriad of additional algorithms have been developed for the analysis of GEMs for strain engineering, contextualization of experimental data, and building cell- and tissue-type specific GEMs [11-13]. In addition to optimization problems that can be solved using linear programming such as FBA, problems have been formulated to take advantage of mixed integer linear programming (MILP; see [14] for a review of optimization problems in systems biology). MILP employs binary state variables during optimization to solve problems that involve discrete activation or inactivation of variables. MILP problems are particularly well-suited to network-based models, since they allow switch-like behavior that can include or exclude network components (e.g., shutting reactions off/on). MILP has been used widely for gap-filling of GEMs, a process in which constraints or objectives are set to recapitulate a known phenotype by adding biochemical functions from a universal set of reactions [15]. In MILP problems used for gap-filling, the objective function is generally minimization of the number of modifications to a GEM that must be made to satisfy the constraints imposed (e.g., metabolite uptake or secretion, production of biomass). One consequence of this formulation is that alternative solutions, which contain unique sets of reactions which need to be altered in the network or added, are common for large networks that have a large space of potential solutions to draw from (e.g., a large universal set of reactions). These alternative optima in MILP problems are increasingly being considered and leveraged to understand redundancy in solutions and whether or not portions of a solution may be spurious [16-18].

It has been shown that the order in which separate instances of gap-filling are applied to the same network (e.g., gap-filling for growth on individual carbon sources iteratively) strongly influences which reactions are included in the resulting network [19]. In this same study, the alternative solutions generated during this process were used to improve gene essentiality predictions using EnsembleFBA, a technique in which sets of alternative GEMs are used to perform FBA to determine gene essentiality. Using the entire ensemble, performance can be tuned by varying the voting threshold required to make a specific prediction. This is analogous to the threshold-based voting procedure used to construct receiver operating characteristic curves for ensemble-based machine learning models such as random forest [20]. This approach is likely to be highly beneficial for studies of organisms for which little biochemical data are available, which typically have many gaps in their GENRE and thus have many highly-variable alternative gap-filling solutions. Although a nascent approach for studying GENREs, we have built on these observations, and ensemble generation and analysis have been applied in several cases [6,19,21,22].

Here, we present Medusa, a Python package for the generation and analysis of ensembles of GENREs. Medusa provides a framework for compactly representing ensembles of GENREs, avoiding the redundancy of storing many separate models while still being flexible enough to represent variation in any component within a GENRE. Medusa manages ensemble storage and indexing during simulation, allowing users to interact with an entire ensemble in the same way they would interact with an individual GENRE using any constraint-based reconstruction and analysis (COBRA) method. Furthermore, by standardizing the representation of ensembles and their interface with existing COBRA methods, Medusa enables the application of supervised and unsupervised machine learning to gain insight into the influence of varying components within an ensemble of GENREs on the predictions they make. The architecture and functionality of Medusa were designed to make ensemble analyses as accessible and usable as COBRA methods applied to single networks.

Design and implementation

Architecture overview and dependencies

Medusa is built on top of COBRApy, a Python-based package in which many COBRA methods are implemented [23]. Although a dependency-free approach in which ensemble simulation methods are implemented from the ground up could be more efficient, we chose to extend COBRApy to greatly decrease the size and complexity of the codebase and to reduce the domain-specific knowledge required to use Medusa and understand the source code (i.e. decrease the effort for existing COBRApy users and contributors to use Medusa). As such, the architecture of Medusa closely mimics COBRApy .

Architecture of ensembles and simulation output in Medusa.

A) Ensemble functionality is implemented in Medusa through the Ensemble class. The Ensemble class exposes three attributes to the user: features, members, and base_model. Both features and members are container-like objects similar to genes, reactions, and metabolites in COBRApy. Within each container, Medusa objects of class feature and member are stored. The base_model attribute points to a COBRApy model. This base_model contains all of the features present in any member of the ensemble, and is manipulated when generating an ensemble or performing simulations. B) Description and shape of simulation results for common COBRA methods as implemented in Medusa. Each simulation method returns a distribution-equivalent to the single-model simulation result.

At the time of this writing, GENREs are represented within COBRApy using a Model class. The Model class manages the interface between COBRApy and numerical solvers through optlang, a Python package for formulating and solving optimization problems that extends the symbolic mathematics package SymPy [24,25]. GENREs are represented by a Model using additional classes with biological analogs (Metabolite, Reaction, and Gene). Objects belonging to each of these classes are stored within container-like objects (metabolites, reactions, and genes, respectively) that are each an attribute of a Model. Each Metabolite, Reaction, and Gene has attributes which might affect simulations performed using the Model, such as the lower and upper bounds of flux through each Reaction or the gene-protein-reaction relationship for each Reaction, which link them to specific Genes.

In Medusa, ensemble functionality is introduced using three new classes. The first, Feature, describes a GENRE component which has a parameter that varies across an ensemble (e.g., a reaction that is reversible in some ensemble members but irreversible in others). The second, Member, describes individual GENREs within an ensemble and their state for each Feature. The third is the Ensemble class, which references every Feature and Member associated with an ensemble of GENREs, as well as a COBRApy Model, referred to as the “base model”. This base model holds all COBRApy objects that might be associated with any Member in the Ensemble. Within an Ensemble, each Feature references the component within the base model (Metabolite, Reaction, or Gene) for which it encodes alternative parameter values, as well as the attribute within that component that is modified (e.g., the upper bound of flux through a reaction). When a simulation is to be performed using a particular GENRE within an Ensemble, Medusa changes the state of the base model to represent the proper state for the corresponding Member for every Feature. Thus, the Ensemble can represent any number of variants in GENRE structure throughout an ensemble (e.g., reaction presence/absence, reversibility, alternative gene-protein-reaction relationships) and can be used to apply any methods implemented in COBRApy. Furthermore, this implementation has a memory footprint only slightly larger than a single COBRApy Model, and facilitates queuing of simulations for parallel processing.

Medusa is developed partially with test-driven development. Unit tests are implemented using the pytest package (https://docs.pytest.org/en/latest) and are run automatically with each modification to the Medusa github repository via continuous integration with TravisCI (Travis CI, GMBH, Berlin, Germany). Support is provided for Python version 3.4 and later. Tests are run in TravisCI using Python 3.6, with plans to include Python 3.7 and 3.8. No support is provided for Python 2.

Results

Performing ensemble simulations

Currently, users can perform FBA, flux variability analysis (FVA), single gene deletions, and single reaction deletions using an ensemble in Medusa. In each case, the simulations are performed with a single function, which returns a collection of results, where each entry corresponds to the simulation results for a single ensemble member. Users have the option of performing simulations using the entire ensemble, a specific set of Members, or a random fraction of Members. For the currently implemented analyses, the collection of results are returned as a DataFrame from the pandas Python package, where each column corresponds to the entry normally populating the results for a single network (e.g., a reaction ID for FBA/FVA, a gene ID for single gene deletions), and each row corresponds to an ensemble Member (except for FVA, where two columns are required for each ensemble member to describe the minimum and maximum fluxes). See for a schematic describing how the shape of data describing simulation results changes for each simulation method.

Usage and benchmarking of Medusa.

A) Comparison of ensemble flux balance analysis (FBA) simulations across different conditions. A) Ensemble FBA performed on glucose or mannose minimal media using an ensemble of 1000 GEMs for Staphylococcus aureus. This ensemble was generated in [6] by iteratively gap-filling a draft reconstruction to enable biomass production in single-carbon source growth conditions supported with experimental data. Mean for the distribution for either condition shown by vertical line of same color. Predicted flux through biomass is higher on mannose than glucose (Wilcoxon signed-rank test, p < 1E-5). B) Benchmarking of Ensemble FBA. Simulations were performed with the same S. aureus ensemble of 1000 members as in panel A. C) Benchmarking of ensemble generation time with iterative gapfilling (process shown in Fig 3). D Comparison of Ensemble FBA simulations before and after curation using a machine learning-guided approach. Ensemble FBA was performed on complete medium (uptake of -1000 mmol/g*DW*hr allowed for all metabolites) using an ensemble of 1000 GEMs for S. aureus. Mean for the distributions before and after curation shown by vertical line of same color in the last panel in the workflow. Ensemble and machine-learning guided curation identified N-Formimino-L-glutamate iminohydrolase as a driver of variation in simulated flux through biomass. Based on a literature search, this reaction was inactivated in all ensemble members and ensemble FBA was performed again, resulting in the shift in the distribution shown. See documentation for full narrative-style example: https://medusa.readthedocs.io/en/latest/machine_learning.html.

Fig 3

Iterative gap-filling strategy used to generate ensembles implemented in Medusa.

A) Given a list of conditions in which an organism satisfied some objective (e.g., secretion of a specific metabolite, growth) and a draft GENRE for the organism, gap-filling is performed sequentially on each condition. After each gap-fill step on a single condition, reactions in the gap-fill solution are added to the GENRE before starting gap-filling on the next condition. In this schematic, all conditions are used, but Medusa also allows users to randomly subsample a fraction of conditions to generate more variation in the resulting gap-filled GENREs. Here, gap-filling is performed to enable growth in the presence of each individual metabolite (indicated by a green arrow for a single condition during each gap-fill step). B) Medusa iteratively performs gap-filling as shown in panel A, then shuffles the order of conditions to introduce variation in the gap-fill solution. After a user-defined number of cycles through the process shown in panel A, Medusa generates an ensemble containing all unique GENREs that resulted from gap-filling.

Because the Ensemble object implemented in Medusa maintains a COBRApy Model object, users can also perform any custom simulation they would like by 1) manipulating the COBRApy Model to be suitable for their simulations (e.g., set custom constraints or objectives), 2) setting the state of the model to represent an ensemble member using Medusa functionality, 3) performing their simulation, then 4) iterating through any other ensemble members they would like and performing steps 2–3.

Comparing ensemble simulations

When performing simulations with an ensemble rather than a single model, results shift from single values to distributions. While this explicitly accounts for uncertainty, it also requires that statistical approaches are applied to interpret differences in distributions. In the simplest case, performing ensemble FBA to predict the growth rate on each of two different media conditions for a single bacterial species generates two distributions of predicted growth rates . A paired univariate test (e.g., t-test or a non-parametric equivalent) can be used to determine whether the predicted growth rate is equal in these two conditions. This example is demonstrated in the Medusa documentation at https://medusa.readthedocs.io/en/latest/stats_compare.html.

Coupling ensemble modeling with machine learning

The availability of ensemble generation and simulation methods in Medusa provides ample opportunity to apply machine learning to leverage variation in ensembles. One application area that we have developed focuses on guiding the curation of genome-scale metabolic network reconstructions by attributing simulation uncertainty to network components in the ensemble, then prioritizing curation of these components based on how much they contribute to simulation uncertainty [6]. This approach can be broken down into four steps: 1) ensemble generation, 2) ensemble simulations, 3) unsupervised learning to summarize simulation uncertainty, and 4) supervised learning to associate uncertainty in network structure with uncertainty in simulations. In addition to our published work utilizing Medusa for this purpose, we provide an example in the Medusa documentation that applies this method to a single ensemble FBA simulation: https://medusa.readthedocs.io/en/latest/machine_learning.html.

In this example, a previously generated ensemble is loaded into Medusa, media conditions are set to allow uptake of any metabolites with transporters in the ensemble, and ensemble FBA is performed. Then, a random forest regressor is used to predict the simulated values for flux through biomass (e.g., the values generated with ensemble FBA) for each ensemble member using the Medusa states (i.e., binary reaction presence/absence) for the same ensemble member as input. Examining the most predictive feature in the random forest regressor, we perform a literature search and find that the feature is likely not present for the bacterial species under study (i.e. the species does not have the ability to catalyze the reaction described by the feature). Based on this examination, we disable the feature (which is a reaction in this case), perform ensemble FBA again, and find that this curation step has reduced the average predicted flux through biomass . Although we use this approach to guide curation, we also envision the same process having great utility for attributing simulation uncertainty more generically, such as with simulations performed using an ensemble generated using ‘omics integration methods [11].

Generating ensembles

The simplest way to generate an ensemble in Medusa is to load a collection of models from files or a collection of COBRApy Model objects. This can be performed by specifying a batch of files that Medusa loads iteratively to construct an ensemble or by loading all the models into memory as COBRApy Models and passing them to an ensemble constructor function. Access to both of these methods for constructing an ensemble gives users flexibility to use other software that generates multiple GENREs and to construct their own ensembles from COBRApy Models that are only ever present in memory (i.e., never written to disk). There are also performance tradeoffs between each method; constructing an ensemble from a batch of files allows the user to specify a batch size for the number of models loaded into memory during each iteration. Increasing the batch size increases memory usage but generally decreases the time required to construct an ensemble. Similarly, constructing an ensemble from a collection of COBRApy Models already loaded into memory has high memory usage but low runtime.

Ensembles can also be generated in Medusa by performing gap-filling on an individual GENRE. Medusa implements a previously-developed algorithm for gap-filling GENREs using growth phenotyping data [6,19]. The algorithm takes a GENRE with an objective function (i.e., a COBRApy Model), a universal reaction database (stored as a COBRApy Model) and a dataset of binary growth/no-growth calls on defined media conditions as input. The objective function is then set as a constraint with bounds such that any feasible flux distribution must enable activity within the bounds (e.g. at least some amount of flux through biomass production or a demand reaction). All reactions from the universal reaction database are added to the GENRE being gap-filled, and a new objective function is set to minimize the sum of fluxes through reactions in the reaction database. The problem is solved to identify reactions taking part in this minimal flux activity. To generate a single gap-filled ensemble member, the draft GENRE is iteratively gap-filled on each positive growth media condition. This process is repeated to produce the number of desired ensemble members. Variation in model structure is introduced by randomizing the order in which media conditions are used for gap-filling. We previously found that ensemble members generated with this iterative strategy are equivalent to gap-filled models generated by finding a single gap-filling solution to all growth conditions simultaneously (i.e., with a single optimization step), yet this iterative process is multiple orders of magnitude faster than the global solution [19]. See for a schematic summarizing this approach.

This gap-filling process is implemented in Medusa through a single function, and a full example of preparing a model and all data necessary for this process are provided in the Medusa documentation at https://medusa.readthedocs.io/en/latest/creating_ensemble.html. Although Medusa implements the previously published version of this approach, it also allows users to randomly subsample a fraction of conditions to generate more variation in the resulting gap-filled GENREs. We benchmarked construction of ensembles using the standard Medusa iterative gap-filling approach, a Staphylococcus aureus draft metabolic network reconstruction generated via ModelSEED [26], and the modelSEED universal biochemistry using a 2019 Apple Macbook Pro (2.4GHz Intel Core i5-8279U, 16GB RAM). Ensembles with 5, 25, or 100 members, gap-filled on 10 randomized positive growth media conditions, took an average of 1.8, 5.3, and 18.6 minutes to generate in 10 independent trials (). A substantial portion of this time (~60 seconds) is spent copying the universal model; time spent after this copy step scales linearly with the size of the ensemble being generated.

Benchmarking Medusa

The structure of a Medusa Ensemble was designed to help users simplify their code when dealing with an ensemble of models, decrease memory and disk storage demands, and decrease processing time when possible. The Medusa documentation provides thorough benchmarking examples for multiprocessing using Medusa, constructing an ensemble, and memory and disk utilization. These benchmarking notebooks demonstrate the memory savings when using a Medusa Ensemble rather than individual models (~50MB with Medusa vs. ~16GB for individual models representing an ensemble with 1000 members) and the disk space savings in the same scenario (~6MB for a Medusa Ensemble vs. 1GB for 1000 individual models). Medusa reduces memory and storage demands by over two orders of magnitude in these cases. Improvements in runtime for performing FBA in Medusa, rather than using individual models, are more modest but improve further when using more than one processor. shows the decrease in optimization runtime in Medusa as the number of processors is increased, and the diminishing returns associated with a larger number of processors. Memory, storage, and runtime benchmarking were performed using the same hardware as the iterative gapfilling benchmarking in . All benchmarking is available in narrative-style notebooks in the Medusa documentation.

Availability and future directions

Stable releases of Medusa are available through the Python package index (PyPI, https://pypi.org) as well as github (https://github.com/opencobra/Medusa). Documentation is available through Readthedocs at https://medusa.readthedocs.io. Current development efforts are focused on parallelization and integrating Medusa with other Python-based tools in the COBRA community.

Currently, the only high-level function available to users of Medusa to share ensembles generates a serialized version of the Ensemble object using the pickle package (the standard Python library package for serializing objects). Alternatively, users can save the base_model for any Ensemble as a Systems Biology Markup Language (SBML, [27]) file using COBRApy, then choose the formatting option of their liking to save Feature and Member information for the Ensemble. In SBML, there is not currently a standardized way to represent ensembles as required in Medusa. We plan on extending the Flux Balance Constraints package, an extension to SBML intended for constraint-based models, to enable standardized sharing of Medusa ensembles. Until then, we recommend users include both an SBML file for the base_model of each Ensemble and the serialized pickle at time of publication.

22 Oct 2019

Dear Jason,

Thank you very much for submitting your manuscript, 'Medusa: software to build and analyze ensembles of genome-scale metabolic network reconstructions', to PLOS Computational Biology. As with all papers submitted to the journal, yours was fully evaluated by the PLOS Computational Biology editorial team, and in this case, by independent peer reviewers. The reviewers appreciated the attention to an important topic but identified some aspects of the manuscript that should be improved.

We would therefore like to ask you to modify the manuscript according to the review recommendations before we can consider your manuscript for acceptance. Your revisions should address all the specific points made by each reviewer and we encourage you to respond to particular issues. In particular, please be sure to address the points about quantification of the performances and gains of MEDUSA.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.raised.

In addition, when you are ready to resubmit, please be prepared to provide the following:

(1) A detailed list of your responses to the review comments and the changes you have made in the manuscript. We require a file of this nature before your manuscript is passed back to the editors.

(2) A copy of your manuscript with the changes highlighted (encouraged). We encourage authors, if possible to show clearly where changes have been made to their manuscript e.g. by highlighting text.

(3) A striking still image to accompany your article (optional). If the image is judged to be suitable by the editors, it may be featured on our website and might be chosen as the issue image for that month. These square, high-quality images should be accompanied by a short caption. Please note as well that there should be no copyright restrictions on the use of the image, so that it can be published under the Open-Access license and be subject only to appropriate attribution.

Before you resubmit your manuscript, please consult our Submission Checklist to ensure your manuscript is formatted correctly for PLOS Computational Biology: http://www.ploscompbiol.org/static/checklist.action. Some key points to remember are:

- Figures uploaded separately as TIFF or EPS files (if you wish, your figures may remain in your main manuscript file in addition).

- Supporting Information uploaded as separate files, titled 'Dataset', 'Figure', 'Table', 'Text', 'Protocol', 'Audio', or 'Video'.

- Funding information in the 'Financial Disclosure' box in the online system.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com  PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

We hope to receive your revised manuscript within the next 30 days. If you anticipate any delay in its return, we ask that you let us know the expected resubmission date by email at ploscompbiol@plos.org.

If you have any questions or concerns while you make these revisions, please let us know.

Sincerely,

Hugues Berry

Associate Editor

PLOS Computational Biology

Douglas Lauffenburger

Deputy Editor

PLOS Computational Biology

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: This software paper describes the Medusa tool, a software to build and analyse ensembles of genome-scale metabolic network reconstructions. The main interest of Medusa is to enable the simultaneous run and comparison of analyses over families of metabolic network, such as FBA, FVA and gene essentiality. To that goal, a family of metabolic networks is represented in CobraPy Python package using a so-called Ensemble class, which references both a CobrayPy Model (describing the base model) and Features and Members to describes the specificities of the family of metabolic networks. This smart representation allow comparing ensemble of simulations and running pipelines which alternate analyses of families of metabolic networks, then selection of sub-families, followed by novel analysis. An example of such a pipeline is a curation procedure based on the prioritizing of components which contribute the most to simulation uncertainty. Another example is the automatic generation of a family of gap-filled metabolic models according to different media conditions allowing the user to see the impact of the order of the selected media conditions over the gapfilling procedure.

As stated by the author, this representation of families of metabolic networks which relevant classes in CobraPy is a very smart solution to run multiple analyses and workflows and study the variability of metabolic networks while avoiding to store the multiple networks generated along the pipelines. This is particularly useful for automatic analyses involving the generation of multiple models.

However, the paper is lacking some elements to ensure that the Medusa package can be used by the community. They correspond to the four main limitations below that needs to be clarified.

(1) There is no information about the gains enabled by Medusa. What is the maximum number of metabolic networks that an Ensemble can describe ? What is the gain of storage space enabled by the ENSEMBLE representation with respect to a usual representation with individual models ? What are the gains of performance for FBA, FVA, and gene essentiality computations with respect to simulations with individual models ?

(2) Few details are given about the generation of ENSEMBLE. How can an ENSEMBLE be generated from a family of SBML files generated with other tools ? Does the package enable this ? In addition, the author suggest to store share ENSEMBLES by using a SBML file for the base_model and then to serialize the ENSEMBLE. What is the gain if this is the recommended solution ?

(3) Figure 4 is quite misleading. What is the advantage of generating Ensembles of gapfilled GENREs by shuffling the different growth conditions ? May the author show FBA and FVA simulations evidencing that the study of the Ensemble generates gains in the study of the considered metabolic network ?

(4) In the introduction, the authors advocate that ENSEMBLEs are useful to handle alternative solutions generated by simulations, including infinite number of feasible flux values within the range determined by a FBA solution space. I do not understand how MEDUSA is able to handle this issue ? Is this possible to generate an ENSEMBLE modeling all such GENREs with alternative feasible fluxes and then to study them with MEDUSA ? How can this be done ? If this is not the case, I would rather suggest to remove the paragraph in the introduction to avoid confusion.

Reviewer #2: The work presented here represents a useful tool for the community. However, there are some points in the manuscript that require clarification.

• In the example provided in the “Coupling ensemble modeling with machine learning,” it is not clear what data the random forest is trained on. Where it says that the “… random forest regressor is used to predict the simulated values …,” are the authors actually training the random forest on those simulations and then going in to the decision trees to determine the most predictive features? If so, this should be clarified. If not, then the information needs to be added.

• Some sense of the simulation “wall time” would be helpful for various size models/ensembles. Of course, it will also be necessary to provide information regarding hardware used in order for the wall time to have any significance.

• One assumes that Medusa is based on Python 3, and it will become obvious once an individual uses the package. However, it would still be worthwhile to mention the version of Python against which this package is built against.

Reviewer #3: In this paper, the authors describe and demonstrate Medusa, a (Python) tool for the analysis of ensemble of metabolic networks. While constraint-based modelling methods have gained immense popularity over the last two decades, they are still beset with a key problem of "degeneracy" in solutions. Furthermore, there exists uncertainty in network structure, in terms of the topology, and many gap-filling algorithms exist to fill gaps (identify missing reactions) in metabolic networks. Building on COBRAPy, Medusa can perform ensembleFBA, to address these and other important issues. Medusa also lends itself to applying machine learning tools for improving curation, and the documentation lists an example of this too.

I found the paper to be very well-written and crisp. It gives a very clear overview of the use of Medusa, and the underlying methods.

Minor Revisions:

1. It would be nice to describe an example of how parallelisation affords advantages, and also a brief discussion of how it scales. It would be advantageous to the users to be well aware of the parallel features of Medusa, which will enable them to work with larger numbers of simulations (or larger networks).

2. While Figure 3 is very good and explanative, I would find it interesting to look at a scatterplot of difference in biomass growth rates vs. similarity of the reaction networks in the ensemble. That is, for every pair of networks in the ensemble, plot a point corresponding to (d(network1, network2), v_bio1-vbio2). Is there a clear correlation? What would the outliers, if any, mean? For distance, a Jaccard distance could be used.

3. An example (both in the manuscript) and documentation, as to how ensemble performance can be studied by plotting ROC curves in Medusa could be described.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

27 Jan 2020

Submitted filename: medusa_reviews.docx

Click here for additional data file.

17 Mar 2020

Dear Jason,

Thank you very much for submitting your manuscript "Medusa: software to build and analyze ensembles of genome-scale metabolic network reconstructions" for consideration at PLOS Computational Biology. Based on the reviews, we are very likely to ultimately accept this manuscript for publication; however, toward that end we ask that you consider a minor revision of the manuscript to address the final recommendation made by reviewer #1 (about the citation of Biggs & Papin, PLOS Comp Bio 2017, see below).

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Hugues Berry

Associate Editor

PLOS Computational Biology

Douglas Lauffenburger

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors have addressed most of my comments. However, they decided to not directly address my comment on the former figure 4 (now figure 3) about the underlying motivations for the choice of the gap-filling method implemented in medusa to generate ensembl. Although I understand their position, I would suggest as a last minor revision to clearly refer to the paper Biggs & Papin, PLOS Comp Bio 2017 when they introduce Figure 3 (paragraph that starts with "Ensembles can also be generated in Medusa by performing gap-filling on an individual GENRE..." to include a summary of the detailed motivation that they provided in the response to reviewed.

Reviewer #2: The authors have addressed all of my concerns.

Reviewer #3: All my comments have been addressed satisfactorily, more or less. I am happy to recommend the publication of the manuscript.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see

17 Mar 2020

Submitted filename: medusa_reviews.docx

Click here for additional data file.

3 Apr 2020

Dear Jason,

We are pleased to inform you that your manuscript 'Medusa: software to build and analyze ensembles of genome-scale metabolic network reconstructions' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Hugues Berry

Associate Editor

PLOS Computational Biology

Douglas Lauffenburger

Deputy Editor

PLOS Computational Biology

***********************************************************

21 Apr 2020

PCOMPBIOL-D-19-01338R2

Medusa: software to build and analyze ensembles of genome-scale metabolic network reconstructions

Dear Dr Papin,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Laura Mallard

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol