- Research Article
- Open Access
MitoCore: a curated constraint-based model for simulating human central metabolism
© The Author(s). 2017
- Received: 2 June 2017
- Accepted: 13 November 2017
- Published: 25 November 2017
The complexity of metabolic networks can make the origin and impact of changes in central metabolism occurring during diseases difficult to understand. Computer simulations can help unravel this complexity, and progress has advanced in genome-scale metabolic models. However, many models produce unrealistic results when challenged to simulate abnormal metabolism as they include incorrect specification and localisation of reactions and transport steps, incorrect reaction parameters, and confounding of prosthetic groups and free metabolites in reactions. Other common drawbacks are due to their scale, making them difficult to parameterise and simulation results hard to interpret. Therefore, it remains important to develop smaller, manually curated models.
We present MitoCore, a manually curated constraint-based computer model of human metabolism that incorporates the complexity of central metabolism and simulates this metabolism successfully under normal and abnormal physiological conditions, including hypoxia and mitochondrial diseases. MitoCore describes 324 metabolic reactions, 83 transport steps between mitochondrion and cytosol, and 74 metabolite inputs and outputs through the plasma membrane, to produce a model of manageable scale for easy interpretation of results. Its key innovations include a more accurate partitioning of metabolism between cytosol and mitochondrial matrix; better modelling of connecting transport steps; differentiation of prosthetic groups and free co-factors in reactions; and a new representation of the respiratory chain and the proton motive force. MitoCore’s default parameters simulate normal cardiomyocyte metabolism, and to improve usability and allow comparison with other models and types of analysis, its reactions and metabolites have extensive annotation, and cross-reference identifiers from Virtual Metabolic Human database and KEGG. These innovations—including over 100 reactions absent or modified from Recon 2—are necessary to model central metabolism more accurately.
We anticipate MitoCore as a research tool for scientists, from experimentalists looking to interpret their data and test hypotheses, to experienced modellers predicting the consequences of disease or using computationally intensive methods that are infeasible with larger models, as well as a teaching tool for those new to modelling and needing a small, manageable model on which to learn and experiment.
- Constraint-based model
- Metabolic network
- Flux balance analysis
- Central metabolism
- Mitochondrial metabolism
Human central metabolism is a large and complex system under sensitive homeostatic control, and its disturbance is causative or associated with many diseases and responses to toxins. However, it is often difficult to relate more than a handful of these changes to their underlying origin or their down-stream impact, due to the highly connected nature of the reactions of central metabolism. Computer models are widely accepted in many fields as a tool to incorporate complexity and simulate changes, allowing predictions to be made and providing a unifying framework to interpret empirical data, especially from large, noisy and incomplete data sets. Yet modelling is treated with scepticism by many biomedical researchers despite their potential broad utility . Simple models of enzyme kinetics (using the assumptions of Henri-Michaelis-Menten kinetics ) are familiar to biomedical scientists, but are impractical for simulations of central metabolism due to every reaction needing parameterisation, alongside the computational expense of solving the large set of differential equations. However, constraint-based models of metabolism used in conjunction with methods such as flux balance analysis  are particularly useful for simulating metabolic changes in large metabolic networks as they can incorporate flexibility, do not require kinetic parameters and are computationally inexpensive. Many genome-scale constraint-based models [4–10] have covered central metabolism and used successfully to model diseases [11, 12]. But these models do not simulate the realistic production rate of ATP (with the recent exception of Recon 2.2 ), a crucial element of modelling central metabolism. Furthermore, the interpretation of simulation results from thousands of reactions is difficult (especially for new-comers). In addition, attempts to simulate diseases can result in the prediction of physiologically improbable reaction fluxes due to erroneous “short-circuits” and energy-generating cycles . These are caused by several common problems including: incorrect parameters for directionality constraints, the assignment of reactions to the wrong cellular compartments, or inaccurate representations of pathways, enzymes, transport steps, prosthetic groups and metabolites. These errors can introduce unrealistic bypasses and shuttles that appear to compensate for a disease state. For example, proton-coupled mitochondrial transporters running in reverse and thus pumping protons that contribute to ATP generation by the mitochondrial ATP synthase, and the confounding of free co-factors with prosthetic groups, especially the flavin adenine dinucleotide of mitochondrial succinate dehydrogenase (SDH) and electron-transferring flavoprotein (ETF), leading to incorrect electron transport between isolated complexes mediated by bound FAD/FADH. These problems are common in genome-scale models that include an initial auto-generation of the reaction network from databases that can include incomplete or incorrect annotation. These issues are particularly acute for modelling mitochondrial metabolism and metabolite transport, as all the current genome-scale models neglect the electrical gradient component (ΔΨ) of the proton motive force (PMF), and the correct proton cost of making ATP by the mitochondrial ATP synthase in animals . It is also sometimes questionable whether the enormous size and complexity of these genome-scale models benefits simulations where only subsystems of cellular metabolism are of interest, such as central metabolism. Furthermore, the scale of genome-scale models can make some techniques computationally infeasible, such as elementary mode analysis , and the longer runtimes for their simulations and the complexity of the results can hinder exploratory analyses and hypothesis testing.
Many problems of genome-scale models can be avoided by using smaller, curated models validated against data from normal and disease metabolism. A more focussed and carefully defined model allows the user to examine and be confident in each reaction, and more clearly elucidate the behaviour of the system, including any short-comings, as well as interpret the results more easily. These ‘core’ models have been shown to be useful in a range of areas [16, 17] and we previously applied this approach for our iAS253 model of the mitochondrion, which we used to simulate metabolic diseases of the tricarboxylic acid cycle . This model was then used as a basis to simulate other disorders including hypoxia during cardiac ischemia , fumarate hydratase deficiency  and common diseases of the mitochondrial electron transport chain . These simulation results were used to generate detailed mechanistic hypotheses for data interpretation and to design further experiments. However, we recognised that this model could be improved upon by constructing a new model encompassing more of central metabolism and explicitly modelling physiochemical features such as the mitochondrial proton motive force. In particular, ease-of-use could be improved by providing extensive annotation of reactions and their parameters.
Here we present MitoCore, a new constraint-based model of central metabolism that addresses these issues and comprehensively expands upon and refines our previous mitochondrial models. This model has been designed to be easy-to-use, includes extensive annotation, has default parameters to simulate human cardiomyocyte metabolism, and is encoded in the widely used SBML format . We anticipate the model will be of great use to those wishing to interpret empirical data by comparing it to simulations of central metabolism and thus investigate predictive models of disease and toxicology.
Building the MitoCore model of human central metabolism
We designed MitoCore as a constraint-based model of central metabolism with two compartments; one representing the cytosol, outer mitochondrial membrane, inter-membrane space and cytosolic side of the inner mitochondrial membrane, and the other the mitochondrial side of the inner membrane and the mitochondrial matrix. The accuracy and utility of a metabolic model is dependent upon the reactions included and the correct partitioning of these reactions between compartments. Thus to create MitoCore, we built a list of candidate reactions to include by considering human reactions in the KEGG , HumanCyc  and BRENDA  databases that use any metabolites involved in central metabolism, and assigned each reaction to the appropriate cellular compartment(s) by assessing the localisation evidence collated in the MitoMiner database  for its catalysing protein. For reactions catalysed by enzymes with a large amount of evidence for mitochondrial localisation but lacking specific evidence for being in the mitochondrial matrix or matrix side of the inner membrane, we applied the principle of metabolite availability . A summary of this localisation evidence is provided in the mitochondrial evidence section of the supplementary annotation file (Additional file 1, ‘Reaction & Fluxes’ worksheet) and consists of confidence scores from the MitoCarta 2  inventory of genes that encode mitochondrial proteins. Reaction directionality was assigned by taking the consensus from annotation in metabolic databases, estimates of Gibbs free energy [28, 29], and general rules of irreversibility . For each reaction extensive additional annotation was recorded including the original KEGG identifiers, EC number, description, gene mappings (both HUGO gene symbol and Ensembl identifiers), and evidence for the gene’s expression in heart and the protein’s mitochondrial localisation.
We next modelled the remaining transport steps that connect the compartments, and the influx and efflux of metabolites to the cytoplasm. We defined four different transport categories. First were transport steps based on the characterised mitochondrial transport proteins. Many carriers can transport a range of related substrates (although with different affinities) and we modelled each metabolite combination as a separate step including counter exchange and proton-coupling. In some cases proton-coupled transporters were represented by two reactions to model separately the forward and reverse directions, and a proton (plus co-transported PMF metabolite) only used for movement down the proton gradient. This prevented transporters being used in reverse to pump protons artificially (and so transferring PMF) and thereby unrealistically contributing to ATP production by the mitochondrial ATP synthase. The second category was for metabolites whose transporters are unidentified. We modelled these steps as uniporters that are not proton coupled. If the metabolite was charged and moving into the mitochondrial matrix we assumed this would impact ΔΨ and co-transported 0.82 PMF metabolite accordingly. The third category was for metabolites that can diffuse across the inner mitochondrial membrane—including oxygen, carbon dioxide and water—and these were modelled as reversible uniport transport steps. Finally we modelled the insertion of lipids into the mitochondrial inner membrane via the flippases, which was coupled to ATP hydrolysis.
To account for the generation of reactive oxygen species (ROS) from the respiratory chain in MitoCore, we defined their production as 0.001% of the flux through complex I—the primary site of ROS generation from the respiratory chain—and hence also reduce the efficiency of proton-pumping [32, 33]. Subsequent reactions model ROS being converted to water at the expense of NADPH.
Various reactions defined in KEGG do not distinguish between free and prosthetic redox cofactors, e.g. FAD/FADH, which can lead to unrealistic electron transport being mediated between prosthetic redox co-factors. To prevent this in MitoCore, we removed prosthetic co-factors from reactions (such as for complex II and the ETF), and instead directly coupled their reactions, e.g. to the reduction of ubiquinone to ubiquinol.
For convenience in setting up simulations, we created four ‘pseudo’ reactions that summarise aspects of MitoCore’s biological activity and can be used by flux balance analysis as objective functions: ATP hydrolysis (representing cellular ATP demand), and the biosynthesis of heme, lipids and amino acids.
To enable comparison of results from MitoCore to those of genome-scale models such as Recon 2 , we re-used identifiers for metabolites and reactions present in the Virtual Metabolic Human database (https://vmh.uni.lu) for MitoCore where possible. However, it was necessary to create 105 new reactions for MitoCore (Additional file 3) that were either absent from the Virtual Metabolic Human database and Recon 2.2  (such as transport steps or compartment specific versions), inaccurately described (such as specifying prosthetic FAD as a free co-factor), or needed to represent new features (such as the proton motive force). New reaction identifiers were appended with the suffix ‘MitoCore’. Pathways represented in MitoCore include glycolysis, pentose phosphate pathway, TCA cycle, electron transport chain, synthesis and oxidation of fatty acids, ketone body and amino acid degradation and cover all parts of central metabolism involved directly or indirectly with ATP production. Finally, MitoCore was extensively tested to ensure: it contained no erroneous energy-generating cycles; its simulation results were consistent with disorders we investigated previously, such as ischemia and mitochondrial diseases [18–21]; each reaction was capable of having a flux (depending on the constraints placed on the cytosolic boundary transport steps) and so ensure it contained no reaction dead-ends.
To enable others to follow our reasoning for assigning reactions to specific compartments and their directionality, we recorded the provenance for reactions (Additional file 1). For example, the directionality evidence section describes why a constraint has been set; a manually evaluated consensus of information from the KEGG , HumanCyc  and BRENDA  databases, general rules of irreversibility , large ΔG values from eQuilibrator  or estimated using a group contribution method , and information from the literature. In cases for which reaction directionality was unclear, it was kept reversible. To address why a reaction has been included in our cardiomyocyte model, we included a heart expression section consisting of RNAseq and immunochemistry expression levels of genes taken from the Human Protein Atlas (version 14) . The spreadsheet also includes gene mappings, identifiers from Recon 2 and KEGG, mitochondrial localisation evidence (as described above) and baseline reaction fluxes when the objective function was maximum ATP production under normal conditions. For distribution, we encoded MitoCore in SBML  (Additional file 4) and produced a companion annotation Excel spreadsheet (Additional file 1).
Simulating cardiomyocyte metabolism using MitoCore and flux balance analysis
Metabolite degradation and ATP yields with MitoCore and Recon 2.2
To demonstrate how different metabolites are degraded appropriately in the MitoCore model and the effect on ATP production of implementing the proton motive force, we performed a series of simulations using particular ‘fuel’ metabolites in isolation. We simulated the oxidation and degradation of glucose, lactate, hexadecanoic acid, hydroxybutanoate, acetoacetate and 20 different amino acids in separate simulations (Additional file 5). The simulations showed that in each case the degradation routes used were biologically plausible. Unsurprisingly fatty acids were the most energy rich fuels (ATP production of 112 μmol/min/g of dry weight), followed by the amino acids tryptophan (43 μmol/min/g of dry weight), isoleucine (38 μmol/min/g of dry weight), leucine (37 μmol/min/g of dry weight) and phenylalanine (36 μmol/min/g of dry weight) and then glucose (33 μmol/min/g of dry weight).
Comparison of maximum ATP yields from different carbon source ‘fuel’ metabolites
MitoCore ATP Yield
Recon 2.2 ATP Yield
Fatty acid hexadecanoate (C16)
Simulating the mitochondrial disease fumarase deficiency with MitoCore and Recon 2.2
To compare the ability MitoCore and Recon 2.2 to model diseases of central metabolism, we performed preliminary simulations of the mitochondrial disease fumarase deficiency (OMIM:606812). Fumarate deficiency is a mitochondrial disease resulting from a defect in the fumarate hydratase gene (ENSG00000091483) that encodes both cytosolic and mitochondrial isoforms, and that converts fumarate to malate—a key step in the TCA cycle. This can leave patients with almost no residual enzyme activity and causes developmental delay, severe mental retardation, seizures and dysmorphic facial features [34–37]. A key diagnostic marker for the disease is the presence of fumarate in the cerebrospinal fluid.
The effect of proton leak on the electron transport chain in MitoCore
Here we present MitoCore, a curated, constraint-based model of human central metabolism, designed as a predictive model of metabolism in disease and toxicology, and for use by a wide range of researchers. It covers all major pathways involved in central metabolism using 407 reactions and mitochondrial transport steps, and 74 transport steps over the plasma membrane. To increase the metabolic flexibility of MitoCore, we included a large number of reactions that were not assigned to classical metabolic pathways, but could have potentially important roles in supporting central metabolism. MitoCore was parameterised and annotated for cardiomyocyte metabolism, which is useful for many types of analyses as the cardiomyocyte can metabolise a wide range of substrates and has reactions common to many other cell types, as well as representing the metabolism of an organ of utmost importance in human health, disease and toxicology. This allows the simulation results to be generalisable, without having features that are particularly cell specific, such as those found only in hepatocytes. However, we have included reactions that are inactivated in the default cardiac model, but can be activated to represent the metabolic capabilities of other cell types, e.g. gluconeogenesis, ketogenesis, β-alanine synthesis and folate degradation. Thus, the model allows biologically relevant flux distributions to be generated ‘out of the box’ without altering the model, while allowing for easy modification to represent metabolism in other cell types.
MitoCore has several unique features. The first is the more accurate partitioning of metabolism between the mitochondrion and the cytosol by using extensive localisation data and annotation. This partitioning has an important impact on model behaviour as the limited transport steps into the mitochondrial matrix result in dramatic differences in metabolite availability in the matrix compared to the cytosol. Therefore, it is important that reactions are assigned to the correct compartment. We achieved this by manually evaluating the subcellular localisation of each reaction’s catalysing protein by using the mitochondrial localisation evidence in the MitoMiner database .
To model oxidative phosphorylation in MitoCore, we devised a new representation of the PMF and mitochondrial respiratory chain—the second unique feature of MitoCore. MitoCore’s representation of the respiratory chain differs in many key aspects to other metabolic models, due to how it models vectorial protons and accounts for both components of the PMF. MitoCore represents the PMF as a metabolite that is co-transported in steps that transport charged metabolites or protons across the inner mitochondrial membrane, such as the reactions of the respiratory complexes, and in proton-coupled and electrogenic transport steps. This separate modelling of vectorial protons using an additional new PMF metabolite enabled the impact of electrogenic transporters on the PMF to be accounted for in flux balance analysis simulations for the first time, and prevented simulation artefacts where (scalar) protons generated or removed from other parts of metabolism allow unrealistically high ATP production.
The third unique feature of MitoCore is the physiochemical modelling of the transport steps that connect the cytosolic and mitochondrial compartments including their impact on the PMF, and the influx and efflux of metabolites over the plasma membrane. 83 transport steps connect the two compartments, of which 30 are modelled on the known transport mechanisms of characterised transport proteins of the inner mitochondrial membrane, whereas the other 53 represent known transport capabilities of the membrane, such as diffusion of small metabolites. 74 transport steps at the cytosolic boundary represent the import and export of metabolites across the plasma membrane, such as oxygen, carbon dioxide, glucose, fatty acids and amino acids.
MitoCore improves over other models by distinguishing between free and prosthetic redox co-factors. This keeps separate the electrons entering the respiratory chain from different sources, which can otherwise become connected via e.g. a shared flavin adenine nucleotide metabolite, and is particularly relevant under perturbed conditions where the erroneous connection of free and prosthetic flavin adenine nucleotides can cause unrealistic electron transport bypasses to occur. This problem is endemic in large-scale models that auto-generate reaction networks directly from metabolic databases without manual curation of reactions.
Other improvements to MitoCore included modelling protein complexes as one reaction rather than a series of linked reactions, as can often be found in metabolic pathway databases and other models. This is particularly important for simulations of gene knockouts, as a gene deficiency will disrupt the activity of the whole protein complex rather than just one of its subunit’s reactions. Further, to facilitate gene-based analyses (such as knock out studies) we provide reaction-to-gene mappings with both gene symbols and Ensembl identifiers. Finally, to summarise different aspects of central metabolism during simulations with MitoCore, we defined four ‘pseudo’ reactions for biosynthesis of amino acids, lipids and heme that are required by cells for maintenance and growth, and ATP hydrolysis to reflect the energy demand of the cell. These pseudo reactions are designed to be used as objective functions during flux balance analysis and can determine if particular functions of central metabolism have become impaired during simulations of disease.
A weakness when it comes to interpreting many metabolic models is the lack of provenance about their components. For example recording why reactions have been included, why directionality constraints have been set, or the origin of reaction parameters. Thus for MitoCore we created a supplementary annotation spreadsheet (Additional file 1) to record this provenance. The spreadsheet also serves as a useful template to map flux distributions onto, as reaction fluxes can be grouped in an intuitive way directly against useful supplementary information. When combined with the small size of the model, simulations can be generated quickly and then easily interpreted.
To demonstrate that MitoCore produces physiologically realistic results, we simulated cardiomyocyte metabolism using the default parameters with flux balance analysis  and the objective of maximum ATP production. The reaction fluxes showed central metabolism was modelled correctly, with the largest fluxes through the TCA cycle and respiratory chain. Numerous fuel sources (fatty acids, glucose, lactate, ketone bodies and amino acids) were imported, degraded and entered the TCA cycle at several different points (Fig. 1). The combination of different fuel metabolites used in the simulations showed similar trends to experimental ranges in well-perfused heart (although the reported ranges are large, reflecting the wide range of metabolites that the heart can potentially use). For example, experimental measurements of the sources of acetyl-CoA for the TCA cycle  report the majority of acetyl-CoA (60–90%) derives from fatty acids, which is close to the 55% in the MitoCore simulation, whereas glycolysis (including lactate oxidation) accounts for 10–40% of acetyl-CoA compared to 35% in the MitoCore simulation, with the remainder of acetyl-CoA derived from amino acids. Finally, to show the robustness of the simulation results, we performed Flux Variability Analysis (FVA). The ranges for reaction fluxes calculated (Additional file 1) showed a realistic range of values and a lack of flux loops (apart from interrelated mitochondrial transport steps), with most of the fluxes uniquely determined or having a limited range with large non-fixed fluxes corresponding to isoenzymes or transport steps whose fluxes exchange while maintaining a fixed total flux.
To show the importance of explicitly modelling the PMF, we simulated the maximum ATP production achievable using 1 μmol/min/g of dry weight of common metabolic fuels in isolation (effectively calculating their ATP yields) and compared the results to those generated by the genome-scale model Recon 2.2  and from theoretical calculations  (Table 1 and Additional file 5). Recon 2.2 is an update of the widely used Recon 2 model , which also separates protons used by the electron transport chain and the aspartate-glutamate carrier from the rest of metabolism to allow it to predict ATP yields, unlike its predecessor, Recon 2, that produces infinite ATP under these conditions . However, Recon 2.2 does not directly capture the energetic cost on ATP production of other proton coupled mitochondrial transport steps or the ΔΨ part of the PMF. In MitoCore each glucose produced 33 ATP in comparison to 32 calculated theoretically  and 32 in Recon 2.2 . For the fatty acid hexadecanonic acid 112 ATP were produced compared to 108 theoretically  and 107 in Recon 2.2. A significant difference between the models is the number of protons required by mitochondrial ATP synthase to produce one molecule of ATP: MitoCore uses 2.7 (based on the structure of the bovine mitochondrial ATP synthase ) whereas the theoretical calculations use 2.5 and Recon 2.2 uses 4.0 (presumably to account indirectly for the electroneutral proton-coupled phosphate carrier and electrogenic charge-coupled ATP/ADP exchange). Due to also considering additional factors that affect ATP production—including the impact and bioenergetic cost of all the transport steps on the PMF as well as ROS production and removal—we believe our figure is likely to be more accurate than both Recon 2.2 and these theoretical calculations. When comparing the flux distributions between MitoCore and Recon 2.2 it appears that despite having similar ATP yields for many metabolites, the pathways used by MitoCore are more biologically reasonable for some metabolites (i.e. using canonical degradation pathways) and lacked some of the unlikely elements found in Recon 2.2, such as using a lactate shuttle instead of the pyruvate carrier during glucose oxidation. In some Recon 2.2 simulations the malate-aspartate shuttle was inactive, presumably due to its use having a direct penalty on ATP production, the only transport step defined in Recon 2.2 to do so. MitoCore avoids this problem of unrealistic bypasses by using PMF penalties on all the relevant mitochondrial transport steps.
Further differences could be seen when comparing the results of preliminary simulations of fumarase deficiency (Fig. 2) that showed MitoCore replicating the efflux of fumarate seen in patients while also using plausible compensatory mechanisms, which would provide a promising starting point for further simulations and analysis (such as determining the effect of varying parameters). Recon 2.2 did not replicate the efflux of fumarate, and ATP production was maintained by using a reaction (CITL, a non-ATP dependent citrate lyase) that instead effluxed acetate. This reaction is unlikely to be present in human, which suggests the Recon 2.2 model would require further investigation and changes to simulate fumarase deficiency realistically.
An important aspect of modelling is exploring different hypotheses by running multiple simulations with different parameters. This is especially true when simulating mechanisms of disease, which may reveal compensatory mechanisms and the effect of supplying different metabolites to compensate for metabolite deficiencies. Therefore, it is beneficial if simulation runtimes are comparatively short. Even when using computationally expensive types of flux balance analysis, such as geometric FBA, MitoCore takes ~10 s per run, compared to ~3 h per run for Recon 2.2. The short runtime for MitoCore simulations makes the exploration of varying parameters much easier and routine, computationally expensive methods feasible, and allows rapid testing of new ideas and hypotheses.
When comparing maximum ATP production using default parameters to the previous iAS253 model using the same parameters , ATP production was notably lower (101 vs 140 μmol/min/g of dry weight) due to the introduction of PMF metabolites in the electron transport chain and transport steps, demonstrating it has an important effect. The results from the proton leak simulations (Additional file 7) show the use of metabolites to represent the PMF enables the model to replicate experimental observations under perturbed conditions that would otherwise be impossible, such as the reversal of mitochondrial ATP synthase, the ATP/ADP carrier and the phosphate carrier. Such behaviour is well-known in the absence of respiratory chain activity, for instance in ρ0 cells to maintain a proton motive force in mitochondria [42, 43].
Taken together these simulation results show the MitoCore model is capable of producing realistic results using a wide range of metabolites while avoiding some of the common problems of genome-scale models, including: unlikely shuttling between compartments, reactions that are not compartmentalized correctly, incorrect directionality constraints and the inclusion of reactions that are unlikely to be present in human. These problems, while understandable due to the scale of these models, nevertheless make producing biological relevant flux distributions difficult, especially under perturbed conditions, and show the continued need for carefully curated ‘core’ models for some types of simulation. Further, these comparisons illustrate an additional role for ‘core’ models, which is the ability to identify and diagnose problem areas for improvement in genome-scale models.
MitoCore is a constraint-based model of human central metabolism, provided with default parameters that provide physiologically realistic reaction fluxes for cardiomyocytes in FBA simulations. The model has several innovations including a new representation of the respiratory chain and proton motive force, and partitioning of reactions to subcellular compartments based on the latest localisation evidence. To achieve an accurate depiction of central metabolism, each of MitoCore’s 407 reactions and transport steps was manually evaluated for directionality, expression of its gene in heart, and subcellular localisation of its protein. To allow MitoCore to be easily used and to make the results comparable with other models and compatible with other types of analyses, we used identifiers from the Virtual Metabolic Human database for both reactions and metabolites where possible, and recorded KEGG identifiers in the annotation. To help ease of use, we also provide an annotation spreadsheet that provides gene mappings, localisation and heart expression evidence, and notes on parameter choice. MitoCore is provided in SBML format to be compatible with a wide range of software. We hope MitoCore will be of use as a research tool to a wide range of biomedical scientists and students—from experienced modellers interested in central metabolism or using computationally intensive methods that are infeasible on a genome-scale model, to those new to modelling who would like to begin by using a small manageable model, with application as a predictive model of disease.
Identifying reactions of central metabolism to include in MitoCore
The starting point for the MitoCore model was an updated version of the iAS253 mitochondrial model . The model was expanded by searching for human reactions in KEGG , HumanCyc  and BRENDA  that were missing from the iAS253 model and could impact central metabolism. Each reaction was reassessed for subcellular localisation and directionality (see below), and to ensure the reactions were occurring in most tissues including cardiac, tissue expression of their genes and proteins was verified by using the Human Protein Atlas .
Partitioning reactions between the cytosol and mitochondrion in MitoCore
To partition the reactions into either the cytosol or mitochondrion, each enzyme was manually evaluated by using the mitochondrial localisation evidence in the MitoMiner database . MitoMiner collates GFP tagging, large-scale mass-spectrometry mitochondrial proteomics studies and mitochondrial targeting sequence predictions, with detailed annotation from the Gene Ontology and metabolic pathway data from KEGG. MitoMiner also contains homology information allowing localisation evidence to be shared amongst species. All available experimental localisation evidence, mitochondrial targeting sequence predictions and annotation were considered, including from homologs from mouse, rat and yeasts. For reactions with strong evidence for mitochondrial localisation, but where matrix localisation is unclear, the principal of metabolite availability was applied —a reaction can only be present in a compartment if all its substrates are available and its products can be used by reactions within the same compartment . Reactions residing in the mitochondrial matrix or matrix side of the mitochondrial inner membrane were assigned to the mitochondrial compartment. Transport steps were created to connect the cytosol and matrix based on the transport properties of the membrane (active transport, diffusion, etc.). Each reaction was cross-referenced with the Virtual Metabolic Human database (https://vmh.uni.lu) and used its identifiers for reactions and metabolites where possible.
Assigning reaction directionality in MitoCore
Reaction directionality was manually evaluated for each reaction in MitoCore. The KEGG , HumanCyc  and BRENDA  databases were consulted, and general rules of irreversibility were taken into account (such as most reactions consume ATP rather than produce it and carbon dioxide is normally produced not consumed) . The ΔG values for reactions were also considered, both calculated by eQuilibrator  and estimated by using the group contribution method , and large changes noted. If reaction directionality was conflicted or unclear, then the literature was consulted. If support for irreversibility was poor or a consensus could not be found, the reaction was assigned reversible. The information used to make each assessment was recorded in the directionality evidence column of the annotation spreadsheet (Additional file 1). Further refinement of reaction directionality was used to eliminate loops that could produce metabolites such as ATP and NADH for ‘free’, and to prevent the interconversion of NADH and NADPH, unless it had been experimentally verified.
Defining and updating reactions in MitoCore
Modelling of the proton motive force (PMF) was improved across the mitochondrial inner membrane, a pseudo-metabolite was introduced to model the effect on the PMF of proton and electrogenic transport steps across the inner mitochondrial membrane. In total 43 reactions representing respiratory complexes as well as proton-coupled transport steps were rewritten to use this new species. Next, to prevent unrealistic bypasses between free and complex-bound prosthetic flavin adenine dinucleotides, FADH and FAD+ were removed from all reactions and replaced with ubiquinone and ubiquinol. Finally, to represent better the current understanding of central metabolism pathways, some reactions were rewritten, such as the generation of ROS from complex I, and the stoichiometry of the proton-pumping of the respiratory complexes. In some cases a reaction’s subcellular localisation was changed by defining mitochondrial-specific metabolite species, or new reactions were written using existing metabolite species.
To highlight where there are differences between a MitoCore reaction and the corresponding reaction in the Virtual Metabolic Human database, the MitoCore identifier used the Virtual Metabolic Human database identifier with the suffix ‘MitoCore’. (N.B. Recon 2 also uses The Virtual Metabolic Human database identifiers.)
Testing the MitoCore model
The model was extensively tested: each reaction in the model was set as the objective function in a series of simulations to ensure all reactions were capable of having flux, with directionality constraints reevaluated for both the reaction and other members of the same pathway if this was not the case; erroneous energy-generating cycles were manually identified and removed by running a succession of simulations where ATP production was maximised and a wide range of different ‘fuel’ metabolites provided; that the model used physiologically characterised pathways under normal conditions and used plausible mechanisms under perturbed conditions.
The MitoCore model was encoded in SBML v2.1 and its validity checked with the SBML online validator .
Metabolism was simulated by using flux balance analysis (FBA), which can be summarised as calculating the reaction turnover, or fluxes (flows) of metabolites through a network of biochemical reactions assuming a pseudo steady state . The fluxes through the network are constrained by the stoichiometry and directionality of the reactions, as well as flux capacity and cytosol boundary uptake ranges. Cytosol boundary transport steps model the import and export of metabolites to the cell, but the overall rate of production and consumption of metabolites is assumed to be zero, hence a pseudo steady state. A metabolic objective function is chosen for each simulation, and FBA used to calculate an optimal set of reaction fluxes that maximise this function. For simulations here, maximum ATP production (by maximising flux through the pseudo reaction of ATP hydrolysis) was used as the objective function because energy generation is a major purpose of central metabolism, particularly in cardiomyocytes. To determine the robustness of the model solutions and identify the most important reactions contributing to the optimal solution, flux variability analysis (FVA) was used. FVA calculates the flux ranges of each reaction that still give the same (maximal) value of the objective function as the optimal solution, or a fraction of it. FVA was applied at 100% and 98% of the optimal ATP production.
For simulations of ATP yield, all cytosol boundary uptake fluxes for metabolites that can be degraded to produce ATP were set to zero, oxygen was increased to 50 μmol/min/g of dry weight (so that limited oxygen availability did not affect results), while other cytosol boundary conditions were unaltered. The uptake flux of each metabolite of interest was then increased to 1 μmol/min/g of dry weight. Geometric FBA simulations were then performed with maximum ATP production as the objective function. (When Recon 2.2 simulations did not obtain convergence, the ‘epilson’ and ‘flexRel’ parameters were relaxed from the defaults.)
For simulations of fumarase deficiency, the reactions corresponding to mitochondrial and cytosolic fumarate hydratase were inactivated in each model, whilst the uptakes values were left at default values. Simulations were performed with geometric FBA and an objective function of maximum ATP production. To enable geometric FBA of the Recon 2.2 model to converge to a solution, the efflux was prevented of acetone (Recon 2.2 reaction ID: EX_acetone), tetradecanoate (Recon 2.2 reaction ID: EX_ttdca), 3-(4-hydroxyphenyl)pyruvate (Recon 2.2 reaction ID: EX_34hpp_), 4-methyl-2-oxopentanoate (Recon 2.2 reaction ID: EX_34hpp_) and bicarbonate (Recon 2.2 reaction ID: EX_hco3).
For simulations of the effect of proton leak on maximum ATP production, the lower bound of the reaction representing the gene UCP2 (MitoCore reaction ID: HtmB_MitoCore) was increased over a series of simulations, thus forcing a minimum flux–representing proton leak—through the reaction.
We wish to thank Lukasz Zielinski and Alexander Smith for their input and testing of the model, and Edmund Kunji for discussions on the transport of metabolites across the mitochondrial membrane.
ACS, FE and AJR were supported by the Medical Research Council, UK. JPM was supported by the Plan cancer 2014–2019 No BIO 2014 06 and the French Association against Myopathies.
Availability of data and materials
All data generated or analysed during this study are included in this published article and its supplementary information files. In addition, the model and annotation file will be available at the MRC Mitochondrial Biology Unit website (http://www.mrc-mbu.cam.ac.uk/mitocore/).
ACS and AJR conceived the study. ACS and JPM devised the methodology. ACS created the model. ACS, FE, JPM, and AJR participated in the curation and testing of the model. All authors contributed to, read and approved the final manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Kirk PDW, Babtie AC, Stumpf MPH. Systems biology (un)certainties. Science. 2015;350:386–8.View ArticlePubMedGoogle Scholar
- Cornish-Bowden A, Mazat J-P, Nicolas S. Victor Henri: 111 years of his equation. Biochimie. 2014;107 Pt B:161–6.View ArticlePubMedGoogle Scholar
- Orth JD, Thiele I, Palsson BO. What is flux balance analysis? Nat Biotechnol. 2010;28:245–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Duarte NC, Becker SA, Jamshidi N, Thiele I, Mo ML, Vo TD, et al. Global reconstruction of the human metabolic network based on genomic and bibliomic data. Proc Natl Acad Sci U S A. 2007;104:1777–82.View ArticlePubMedPubMed CentralGoogle Scholar
- Thiele I, Swainston N, Fleming RMT, Hoppe A, Sahoo S, Aurich MK, et al. A community-driven global reconstruction of human metabolism. Nat Biotechnol. 2013;31:419–25.View ArticlePubMedGoogle Scholar
- Ma H, Sorokin A, Mazein A, Selkov A, Selkov E, Demin O, et al. The Edinburgh human metabolic network reconstruction and its functional analysis. Mol Syst Biol. 2007;3:135.View ArticlePubMedPubMed CentralGoogle Scholar
- Mardinoglu A, Agren R, Kampf C, Asplund A, Nookaew I, Jacobson P, et al. Integration of clinical data with a genome-scale metabolic model of the human adipocyte. Mol Syst Biol. 2013;9:649.View ArticlePubMedPubMed CentralGoogle Scholar
- Mardinoglu A, Agren R, Kampf C, Asplund A, Uhlén M, Nielsen J. Genome-scale metabolic modelling of hepatocytes reveals serine deficiency in patients with non-alcoholic fatty liver disease. Nat Comms. 2014;5:3083.View ArticleGoogle Scholar
- Uhlén M, Fagerberg L, Hallström BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Tissue-based map of the human proteome. Science. 2015;347:1260419–9.Google Scholar
- Swainston N, Smallbone K, Hefzi H, Dobson PD, Brewer J, Hanscho M, et al. Recon 2.2: from reconstruction to model of human metabolism. Metabolomics. 2016;12:109.View ArticlePubMedPubMed CentralGoogle Scholar
- Yizhak K, Chaneton B, Gottlieb E, Ruppin E. Modeling cancer metabolism on a genome scale. Mol Syst Biol. 2015;11:817.View ArticlePubMedPubMed CentralGoogle Scholar
- Bordbar A, Monk JM, King ZA, Palsson BO. Constraint-based models predict metabolic and associated cellular functions. Nat Rev Genet. 2014;15:107–20.View ArticlePubMedGoogle Scholar
- Fritzemeier CJ, Hartleb D, Szappanos B, Papp B, Lercher MJ. Erroneous energy-generating cycles in published genome scale metabolic networks: identification and removal. PLoS Comput Biol. 2017;13:e1005494.View ArticlePubMedPubMed CentralGoogle Scholar
- Watt IN, Montgomery MG, Runswick MJ, Leslie AGW, Walker JE. Bioenergetic cost of making an adenosine triphosphate molecule in animal mitochondria. Proc Natl Acad Sci U S A. 2010;107:16823–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Trinh CT, Wlaschin A, Srienc F. Elementary mode analysis: a useful metabolic pathway analysis tool for characterizing cellular metabolism. Appl Microbiol Biotechnol. 2009;81:813–26.View ArticlePubMedGoogle Scholar
- Zielinski DC, Jamshidi N, Corbett AJ, Bordbar A, Thomas A, Palsson BO. Systems biology analysis of drivers underlying hallmarks of cancer cell metabolism. Sci Rep. 2017;7:41241.View ArticlePubMedPubMed CentralGoogle Scholar
- Di Filippo M, Colombo R, Damiani C, Pescini D, Gaglio D, Vanoni M, et al. Zooming-in on cancer metabolic rewiring with tissue specific constraint-based models. Comput Biol Chem. 2016;62:60–9.View ArticlePubMedGoogle Scholar
- Smith AC, Robinson AJ. A metabolic model of the mitochondrion and its use in modelling diseases of the tricarboxylic acid cycle. BMC Syst Biol. 2011;5:102.View ArticlePubMedPubMed CentralGoogle Scholar
- Chouchani ET, Pell VR, Gaude E, Aksentijević D, Sundier SY, Robb EL, et al. Ischaemic accumulation of succinate controls reperfusion injury through mitochondrial ROS. Nature. 2014;515:431–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Ashrafian H, Czibik G, Bellahcene M, Aksentijević D, Smith AC, Mitchell SJ, et al. Fumarate is cardioprotective via activation of the Nrf2 antioxidant pathway. Cell Metab. 2012;15:361–71.View ArticlePubMedPubMed CentralGoogle Scholar
- Zieliński ŁP, Smith AC, Smith AG, Robinson AJ. Metabolic flexibility of mitochondrial respiratory chain disorders predicted by computer modelling. Mitochondrion. 2016;31:45–55.View ArticlePubMedPubMed CentralGoogle Scholar
- Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, et al. The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics. 2003;19:524–31.View ArticlePubMedGoogle Scholar
- Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2016;44:D457–62.View ArticlePubMedGoogle Scholar
- Romero P, Wagg J, Green ML, Kaiser D, Krummenacker M, Karp PD. Computational prediction of human metabolic pathways from the complete human genome. Genome Biol. 2005;6:R2.View ArticlePubMedGoogle Scholar
- Chang A, Schomburg I, Placzek S, Jeske L, Ulbrich M, Xiao M, et al. BRENDA in 2015: exciting developments in its 25th year of existence. Nucleic Acids Res. 2015;43:D439–46.View ArticlePubMedGoogle Scholar
- Smith AC, Robinson AJ. MitoMiner v3.1, an update on the mitochondrial proteomics database. Nucleic Acids Res. 2016;44:D1258–61.View ArticlePubMedGoogle Scholar
- Calvo SE, Clauser KR, Mootha VK. MitoCarta2.0: an updated inventory of mammalian mitochondrial proteins. Nucleic Acids Res. 2016;44:D1251–7.View ArticlePubMedGoogle Scholar
- Jankowski MD, Henry CS, Broadbelt LJ, Hatzimanikatis V. Group contribution method for thermodynamic analysis of complex metabolic networks. Biophys J. 2008;95:1487–99.View ArticlePubMedPubMed CentralGoogle Scholar
- Flamholz A, Noor E, Bar-Even A, Milo R. eQuilibrator—the biochemical thermodynamics calculator. Nucleic Acids Res. 2012;40:D770–5.View ArticlePubMedGoogle Scholar
- Ma H, Zeng A-P. Reconstruction of metabolic networks from genome data and analysis of their global structure for various organisms. Bioinformatics. 2003;19:270–7.View ArticlePubMedGoogle Scholar
- Mitchell P. Coupling of phosphorylation to electron and hydrogen transfer by a chemi-osmotic type of mechanism. Nature. 1961;191:144–8.View ArticlePubMedGoogle Scholar
- Cochemé HM, Murphy MP. Complex I is the major site of mitochondrial superoxide production by paraquat. J Biol Chem. 2008;283:1786–98.View ArticlePubMedGoogle Scholar
- Murphy MP. How mitochondria produce reactive oxygen species. Biochem J. 2009;417:1–13.View ArticlePubMedGoogle Scholar
- Phillips TM, Gibson JB, Ellison DA. Fumarate hydratase deficiency in monozygotic twins. Pediatr Neurol. 2006;35:150–3.View ArticlePubMedGoogle Scholar
- Whelan DT, Hill RE, McClorry S. Fumaric aciduria: a new organic aciduria, associated with mental retardation and speech impairment. Clin Chim Acta. 1983;132:301–8.View ArticlePubMedGoogle Scholar
- Kerrigan JF, Aleck KA, Tarby TJ, Bird CR, Heidenreich RA. Fumaric aciduria: clinical and imaging features. Ann Neurol. 2000;47:583–8.View ArticlePubMedGoogle Scholar
- Bourgeron T, Chretien D, Poggi-Bach J, Doonan S, Rabier D, Letouzé P, et al. Mutation of the fumarase gene in two siblings with progressive encephalopathy and fumarase deficiency. J Clin Invest. 1994;93:2514–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Strittmatter L, Li Y, Nakatsuka NJ, Calvo SE, Grabarek Z, Mootha VK. CLYBL is a polymorphic human enzyme with malate synthase and β-methylmalate synthase activity. Hum Mol Genet. 2014;23:2313–23.View ArticlePubMedGoogle Scholar
- Laskowski KR, Russell RR. Uncoupling proteins in heart failure. Curr Heart Fail Rep. 2008;5:75–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Stanley WC, Recchia FA, Lopaschuk GD. Myocardial substrate metabolism in the normal and failing heart. Physiol Rev. 2005;85:1093–129.View ArticlePubMedGoogle Scholar
- Nelson DL, Cox MM, Lehninger AL. Lehninger principles of biochemistry. New York: W. H. Freeman; 2013.Google Scholar
- Dupont CH, Mazat JP, Guerin B. The role of adenine nucleotide translocation in the energization of the inner membrane of mitochondria isolated from ϱ+ and ϱo strains of Saccharomyces Cerevisiae. Biochem Biophys Res Commun. 1985;132:1116–23.View ArticlePubMedGoogle Scholar
- Buchet K, Godinot C. Functional F1-ATPase essential in maintaining growth and membrane potential of human mitochondrial DNA-depleted ρ° cells. J Biol Chem. 1998;273:22983–9.View ArticlePubMedGoogle Scholar
- Bornstein BJ, Keating SM, Jouraku A, Hucka M. LibSBML: an API library for SBML. Bioinformatics. 2008;24:880–1.View ArticlePubMedPubMed CentralGoogle Scholar
- Schellenberger J, Que R, Fleming RMT, Thiele I, Orth JD, Feist AM, et al. Quantitative prediction of cellular metabolism with constraint-based models: the COBRA toolbox v2.0. Nat Protoc. 2011;6:1290–307.View ArticlePubMedPubMed CentralGoogle Scholar