BACKGROUND
metabolic reconstruction is the computational-based process that aims to elucidate the network of metabolites interconnected through reactions catalyzed by activities assigned to one or more genes
<cit> .
the reconstruction process begins with the functionally annotated genome of an organism.
then, the identification of those genes whose putative products catalyze some biochemical reaction, i.e.
gene products with an assigned enzyme commission number  or transport commission number , should be done.
this relational information can be organized in the so-called gene-protein-reaction association tables
<cit> .
in a further step gpr tables will be used to infer candidate metabolic pathways coded in the organism’s genome
<cit> .
in order to automate the process of reconstruction of a target organism’s metabolic network, computational methods have been previously developed that will yield a first draft
<cit> .
this first draft can be used to formulate a mathematical representation of an organism’s metabolism, termed as genome-scale model .
the constraint-based-modeling  is an approach that combines the stoichiometric analysis with optimization techniques to study genome-scale models
<cit> .
the cbm has been successfully used to predict metabolic capabilities such as growth rates, as well as systems responses to environmental or genetic perturbations
<cit> .
when applying cbm to an initial draft of a metabolic model, it is usual to find inconsistencies that can have different causes.
in the initial stages of a metabolic reconstruction, due to annotation errors, as well as the existences of unknown enzyme functionality, gpr associations can be incorrectly established.
thus, some reactions may be not included in the model draft.
as a consequence, some metabolic pathways will contain gaps that will create dead-end metabolites
<cit> .
these metabolites appear in the model as only produced or only consumed by reactions, and hence will never reach a steady state different than the trivial, and then they will never participate in a feasible solution.
they will in turn block any reaction in which they are involved.
there are two classes of dead-end metabolites: i) root-non-produced metabolites  i.e.
metabolites that are only consumed by system’s reactions, and ii) root-non-consumed  that includes those metabolites that are only produced by the network but never consumed
<cit> .
detection of rnp and rnc can be conducted by simply scanning the rows of the stoichiometric matrix.
however, the absence of flow through metabolites rnp  could be propagated downstream  by blocking reactions and thus, additional metabolites would become gaps .
those metabolites that become a gap as a consequence of some rnp metabolite are termed downstream-non-produced .
in a symmetric way, upstream-non-consumed  metabolites are defined as those metabolites that became a gap as a consequence of some present rnc metabolite
<cit> .
in general, the detection of dead-end metabolites and blocked reactions is referred commonly as the gap finding problem
<cit> .
model’s gaps can be solved by adding one or more reactions that allow connecting a dead-end metabolite with other metabolites of the network, a process known as gap-filling
<cit> .
in some cases, the incorporated reactions can be mapped into some coding gene.
however, there are some situations where even if a set of candidate reactions that fill the gaps has been successfully predicted, it could not be possible to find the genes that code for these activities.
in such cases the reactions are called orphan reactions.
methods to predict candidate genes to be assigned to orphan reactions have also been previously developed
<cit> .
thus, the reconstruction of metabolic models is an iterative process in which the cbm plays an important role to detect inconsistencies that should be resolved or curated in order to improve model formulation
<cit> .
automated methods for model curation have been previously developed .
in order to solve the gap-filling problem, an optimization-based method to identify the minimum number of reaction to be included in the model has been proposed by different authors
<cit> .
these methods rely in the use of mixed integer linear programming  combined with universal reaction databases such as kegg
<cit> , bigg
<cit>  or metacyc
<cit> .
other proposed approaches are based on the use of experimental information to detect inconsistencies with model predictions that may suggest errors in model formulation
<cit> .
even though automated methods for metabolic network curation are of an undoubted help to improve model formulation, there may be situations where the manual inspections of a curator are still needed.
this is certainly the case of the reconstruction of networks from genomes that suffer reductive evolution  and code for minimized metabolisms.
during the establishment of symbiosis, metabolic redundancies with the host can result in the loss of enzymatic steps in the endosymbiont network, leading to the emergence of obligate metabolic complementation.
these shared metabolic abilities take the form of interrupted pathways when the endosymbiotic network is reconstructed.
thus, the problem of gap-filling is a complex decision making process where a visual representation of the inconsistencies can help to a model’s curator to understand how gap metabolites and blocked reactions are related and thus find the nature of these inconsistencies.
in this paper we present a method that combines the cbm with an algorithm to compute connected components over bipartite-graphs.
the presented method allows the detection of the isolated sets of blocked and gap metabolites and the way in which these are interconnected in what we have termed unconnected modules  .
then, the analysis of each individual unconnected module simplifies and clarifies the visual representation, and hence can be used to decide how gaps should be filled during the curation of the model.
the availability of accurate gsm is especially relevant in the case of bacterial obligate  endosymbionts since it is not possible the culture of these microorganisms and hence there is an intrinsic difficulty for obtaining direct experimental data from the system.
in this case, modeling can serve as an appropriate proxy for functional characterizations.
for instance, metabolic modeling has been successfully used in the case of buchnera aphidicola  to evaluate the role of the endosymbiont and the host in nitrogen metabolism
<cit> , in sodalis glossinidius  to characterize intermediate steps during the reduction of the network as a result of the interaction with the host metabolism
<cit> , and in blattabacterium cuenoti  to better understand the striking conservation of the endosymbiont metabolism along the evolutionary time as well as its putative role in the nitrogen economy of the system
<cit> .
in this paper we have used the gsm from b. cuenoti  to test the proposed method of curation.
our study allows upgrading icg <dig> to a more accurate model version named imp <dig>  and it can be used as a guide for systems biology experimental explorations of the interaction between cockroaches and their endosymbionts.
methods
constraint-based modeling
the study of the structural properties of biochemical reactions network relies on the analysis of the stoichiometric matrix
<cit> .
let denote by n the stoichiometric matrix associated to a certain metabolic network with m rows and n columns corresponding to the number of metabolites and reactions, respectively.
in the following i and j will refer to the set of metabolite indexes  and reaction indexes , respectively.
moreover, the set of reaction indexes j will be partitioned into two disjoint subsets: the set jint which contains the indexes of internal fluxes, i.e.
the biochemical reactions that take place inside the cell, as well the transport reactions that operate between the cell and the surrounding medium.
on the other hand, the set jex contains the indexes of the exchange fluxes, which are the auxiliary variables used to represent the rate at which certain metabolites are consumed or produced by the system
<cit> .
there will be only one exchange flux per metabolite and these fluxes will be associated to the metabolites belonging to the extra-cellular compartment.
by convention, the activity of the exchange fluxes is defined as positive or negative if the metabolite is produced or consumed by the system, respectively
<cit> .
the cbm approach relays on the use of different kinds of constraints represented by mathematical equations to define the so called flux space f, i.e.
the set of all flux distributions compatible with the given constraints
<cit> .
the steady state condition is imposed over the mass balance equation of each metabolite of the network yielding the following homogeneous system of linear equations:
n.v= <dig>
where the vector v is a flux distribution compatible with the steady-state condition.
moreover, lower and upper bounds are imposed over each reaction to represent additional constraints.
for instance, the thermodynamic constraints that make some reactions to be irreversible are represented by setting to zero the lower bound.
besides, the surrounding environment of a metabolic system can be modeled by setting bounds over the exchange fluxes.
for example, if a given metabolite is available in the medium and thus can be consumed by the system, the lower bound of the corresponding exchange flux should have a negative value.
the lower and upper bounds imposed over each flux can be written as the following systems of linear inequalities:
vjlb≤vj≤vjub∀jεj
where vj is the activity through the flux j, whereas
vjlb and
vjub are its lower and upper bounds, respectively.
together, the homogeneous system of linear equations  and the system of linear inequalities  yields the mathematical representation of the flux space f, expressed as:
f=vεrn:n.v= <dig> vjlb≤vj≤vjub∀jεj
blocked reactions
a reaction in a metabolic model is defined as blocked under a given medium condition if it cannot display a steady-state flux other than zero:
jεjblocked⇔vj= <dig> ∀vεf
where jblocked is the set of blocked reaction indexes.
this set can be computed solving a set of linear programs, as proposed by burgard et al.
<cit> .
the approach consists in calculating the minimum and maximum flux value through each reaction of the system.
when the maximum and minimum values found for a given reaction are both equal to zero, the reaction is said to be blocked under the defined medium condition.
the formulation of the set of linear programs is the following:
min/max:vj∀jεjs.t.n⋅v=0vjlb≤vj≤vjub∀jεj
gap metabolites
gap metabolites in gsm are defined as those vertexes of the network through which there can be no steady state flow
<cit> .
while the detection of rnp and rnc metabolites  is straightforward by scanning of each row of the stoichiometric matrix n, the case of detecting unp and dnc metabolites cannot be accomplished by a simple inspection of the entries of n <cit> .
however, based on the definition of blocked reaction  we found a way to define gap metabolites that allow its identification in a straightforward manner.
definition: a metabolite ί ∈ i in a network under steady-state is a gap if and only if all the reactions in which its participate  belongs to the set of blocked reactions j
blocked
.
this statement implies that there cannot exists a stationary flow through this metabolite.
thus, if we name the set of reactions in which a metabolite ί participate as:
σi=jεj:nij≠0∀iεi
then, the set of gap metabolites igap ⊆ i can be defined as follows:
i∈igap⇔σi⊆jblocked
hence, the detection of gap metabolites can be accomplished by finding the set jblocked and applying equations  and  for each metabolite.
the given definition for a gap metabolite doesn’t make any distinction between the different classes of gap defined .
although there is no general procedure for the classification of gap metabolites as rnp, rnc, unp, or dnc once they have been found, this is possible in simple cases.
for example, if a gap metabolite ί is involved only in an irreversible reaction by which it is consumed, or if all the reactions in which this metabolite is involved are irreversible and in all it is consumed, then metabolite will be a rnp.
a similar reasoning can we argued for the symmetric case, leading to the identification of a rnc.
for more complex cases, a visual inspection of the underling "unconnected module"  may help to the classification.
the coenzyme pseudo-gap problem
in general gap metabolites can be identified using its relation to the set of blocked reactions jblocked as it was explained in previous subsection.
however, there could be some special metabolites that are not "gap" under the definition given by .
nevertheless these metabolites may be the cause that certain reactions get blocked, as the example depicted in figure
<dig> shows.
the metabolite d is not a gap but the mass balances equation for d* implies that v <dig> is equal to v <dig>  and the mass balance equation for d implies v <dig> + v <dig> = v <dig>  as a consequence of these relations v <dig> is restricted to zero, i.e.
v <dig> is a blocked reaction.
this effect propagates downstream to v <dig>  v <dig> and v <dig>  this is why we name d as a pseudo-gap metabolite.
in order to unblock v <dig>  which in term will unblock v <dig>  v <dig> and v <dig>  there must be added a sink for d or d*.
<dig>
these kinds of situations may take place when the biosynthetic pathway of a coenzyme is present in a model, but there are no fluxes draining or degrading the coenzyme produced by this pathway.
these metabolites may be involved in conserved moieties, and in such cases they will be consumed and regenerated in a cyclic manner.
as a consequence they could not be detected as gap metabolites because they will participate in at least two active reactions.
however, if the biosynthetic pathway for a coenzyme is included in a metabolic model, the net production of a coenzyme will not occur under steady state unless some flux consumes it.
hence, the reactions involved in the biosynthesis pathway may become blocked.
a common approach to solve this problem is to include the coenzyme-like metabolite into the biomass equation or alternatively to introduce an exchange flux that can drain the metabolite out of the system.
either of these two situations is equivalent to the addition of the above-mentioned sink.
the way to detect the pseudo-gap metabolites may be summarized as follow.
when an um overlap with the biosynthetic pathway of a certain cofactor, two different situations can be found: the cofactor is included in the um as a gap or it is not.
if it is a gap, the cofactor must be included in the biomass equation in order to solve the um .
if the cofactor is not included in the um , then it is quite probable that the cofactor may be involved in a conservation relation which connect the cofactor to active reactions and for this reason is not detected as a gap under the definition given by equation .
however, even not detected as gap, the cofactor could be the underlying cause that the um become blocked, and for this reason was termed the cofactor a "pseudo-gap" as has been pointed out at the beginning of this section .
the conservation relations can be detected by the analysis of the conserved moieties, which were calculated as previously described
<cit> .
finally, by adding the corresponding cofactor to the biomass equation the um can be solved.
identification of unconnected modules
when analyzing the set of blocked reactions and gap metabolites of a metabolic model, it is common to find relations between both sets due to the fact that blocked reactions may be connected with other blocked reactions through gap metabolite.
in some cases, blocked reactions are directly connected to a rnp/rnc metabolite.
however, this may be not the case of other blocked reactions that may be not linked directly to a rnp or rnc metabolite.
for example, as it can be seen in figure
1a, reaction v <dig>  although blocked, is not directly connected to a rnp metabolite , but indirectly through metabolite through a set dnc metabolites  and other blocked reactions .
similarly occurs in the symmetric case .
as a consequence of these relations, it is possible to systematically establish how the blocked reactions are connected through the gap metabolites.
any metabolic network can be represented as a directed bipartite graph by considering the sets of vertex v = iuj <cit> .
then, a directed edge or arc will exist between a metabolite and a reaction if the metabolite participates in the reaction.
the direction of the arc will be incident to the reaction if the metabolite is a reactant and incident to the metabolite if it is a product.
in the following, the graph associated to a metabolic network will be referred as the metabolic graph.
once the metabolic graph is constructed it is possible to consider any possible sub-graph by selecting a pair of subsets i' ⊆ i and j' ⊆ j of metabolites and reactions respectively.
in particular, it is possible to consider the sub-graph defined by the subset of vertex v = igapujblocked.
this sub-graph will contain the relation that exists between gap metabolites and blocked reactions.
moreover, the set of connected components can be computed over this graph.
in this context, each connected component can be interpreted as a "module" of the metabolism that becomes inactive or unconnected, possibly as a consequence of model inconsistencies, such as the presences of a set of rnp/rnc metabolites.
for this reason the set of connected component will be referred as an unconnected modules .
in simple cases, the reason that causes a certain um to be unconnected from the rest of the network may be found by visual inspection of the graph representing the um.
in such cases connectivity restoring of the rnc and rnp metabolites present in the um may, in general, solve the problem of the unc and dnp metabolites.
the set of elementary operations that can be applied to solve ums has been discussed by kumar and collaborators
<cit> , and it includes: addition of biochemical reactions of transport, incorporation of exchange fluxes and relaxation of irreversibility constraint of some reactions.
more complex situations may include cases such as the pseudo-gap problem described in the previous section.
flux balance analysis
the flux balance analysis  is an approach that combines the description of the flux space f defined by equation  with optimization techniques to find a flux distribution v that maximizes the growth rate
<cit> .
this problem can formulated as a linear program and can be solved with standard techniques of linear programming .
max:vbiomasss.t.n⋅v=0vjlb≤vj≤vjub∀jεj
where the flux vbiomass represents the growth rate of the organism.
this flux, also refereed as the biomass equation, includes all the metabolites that are biomass components, in its specific proportions
<cit> .
in-silico knockout experiments
the fragility analysis of a network was performed by simulating knockout experiments for each metabolic gene included in a gsm.
an in-silico knockout experiment for a given gene consists in bound to zero the flux for each reaction coded by the gene, which are inferred through the grp association table.
after so, fba is used to find the maximal value of the biomass reaction under the genetic perturbation.
if the optimal value of biomass reaction is lower than a certain threshold, the knockout is said to be lethal , otherwise the gene is not essential.
the procedure is performed over all the genes in the model.
minimal medium prediction
the minimal medium may be defined as the smaller set of metabolites that should be present in the medium condition in order a feasible flux distribution v to exist, with a biomass production rate greater than zero.
the minimal medium was calculated by solving a milp algorithm as previously described
<cit> .
the milp algorithm is the following:
min:∑jεjexyjs.t.n⋅v=0vjlb≤vj≤vjub∀jεjintvbiomass≥vbiomasslbvjlbyj≤vj∀jεjexyjε <dig> ∀jεjex
the algorithm requires the incorporation of a set of binary yj variables, one for each exchange flux j ε jex.
moreover, a set of constraints that relates each binary variable with its corresponding exchanges flux should be incorporated to the problem.
then, whenever a binary variable takes the zero value, the corresponding exchange flux is constrained to take a value greater or equal than zero.
additionally, the lower bound corresponding to the biomass production flux should be set to a cutoff value greater than zero in order to guarantee a positive growth rate.
finally, the optimization target is defined in such a way that minimizes the number of active exchange fluxes with a negative value.
due to the fact that each exchange flux is related with one extracellular metabolite this is equivalent to find the minimal set of metabolites that the system must consume in order to produce a biomass flux greater than zero.
detection of reaction subsets
a reaction subsets
<cit>  or full coupling sets <cit>  in a stoichiometric network is a group of reactions that operate together in fixed flux proportions for any flux distribution.
due to the fact that the relations between enzymes and reactions are not always biunivocal, the reaction subset does not always match the concept of enzyme subset previously introduce by pfeiffer et al.
<cit> , and for this reason the term reaction subset seams more appropriate
<cit> .
the rss are structural invariant of network and for that they are independent on the kinetic parameters of the system.
moreover, they shed valuable information that may help to understand how the network is regulated
<cit> .for these reasons the concept rs is important for the analysis of metabolic networks.
in order to compute reaction subset the following pre-processing steps were applied to the network.
first, rows and columns corresponding to gap metabolites and blocked reaction respectively were removed from the stoichiometric matrix.
as it was previously described
<cit> , the constant biomass composition imposed by the stoichiometry of the biomass equation was relaxed by removing the corresponding column, while allowing each biomass component to be drained from the system in an independent way.
then, the identification of rs was done by using the algorithm described in
<cit> .
computational tools
constraint-based analysis was performed using the python-based toolbox cobrapy
<cit> .
lp and milp problems were solved using the gurobi solver
<cit>  acceded through cobrapy.
identification of conservation relations was done by computation the set of extreme rays using the polco package
<cit> .
identification of rs was done by using an implementation based on python
<cit>  of the algorithm described in
<cit> .
the detection of the connected components of a graph was done using the implementation of the algorithm available in the igraph library
<cit> .
graphs were drawn using the yed graph editor
<cit> .
all computation was done on a desktop computer with an intel® core™ i <dig> cpu  <dig> processor, with  <dig>  gib, running under fedora  <dig> linux os.
RESULTS
the first step in the analysis of the gsm of b. cuenoti icg <dig> was to find the sets jblocked and igap of blocked reactions and gap metabolites, respectively.
the results showed that  <dig> reactions over a total of  <dig>  are blocked under any medium condition.
using this information a set of  <dig> metabolites over a total of  <dig>  were detected as gaps.
a bipartite graph representation of the metabolic network was constructed, and the sub-graph defined by the subset of vertex igap u jblocked was selected .
computation of connected components over this sub-graph allows identifying  <dig> um.
after so, each gap metabolite was classified in one of the four categories: rnp, rnc, unc and dnp.
a description of each different um found in icg <dig> is summarized in table
<dig>  which was sorted according to the number of reactions included in each um.
um identified in b. cuenoti icg <dig> gsm.
metabolite abbreviations: mev , 2ombzl , 2ommbl , hxan , xan , r1p , 2dr1p , ura  1pyr5c , 4hthr , u3hga -d-glucosamine), uppg <dig> , apoacp .
the relation between an um and a subsystem was established according to the most frequent subsystem associated to the reactions of the um.
when analyzing the reactions participation of each um, it was found that in many of the cases the set of reactions belonging to an um overlapped with known biochemical subsystems .
for example, all the reaction included in um <dig> belongs to the menaquinol biosynthetic pathway.
it is worth to note that in some cases ums could be composed by an isolated reaction as the cases of um <dig> and um <dig>  after the identification of all um applying the proposed approach, each sub-graph was drawn independently.
visual inspection of each graph was done to detect possible inconsistencies in the network formulation that may help to carry out the manual curation of a gsm.
um <dig> - menaquinol biosynthetic pathway
the biggest um found in the metabolic model  was first analyzed.
figure
<dig> shows the graph that represents the um.
a visual inspection of these figure shows that um <dig> includes all the reactions and metabolite that conforms the menaquinol biosynthetic pathway.
the figure also shows that from the  <dig> gap metabolites included in the um, two are rnp metabolites , and one is an rnc metabolite , while the remaining metabolites are then upstream/downstream gap metabolites.
mevalonate is a precursor for the biosynthesis of isopentenyl diphosphate that in turn is involved in the biosynthetic pathway of many cofactors such as menaquinol and 2-demetyl-menaquinol.
an inspection of the current genome annotation of b. cuenoti <cit>  was done to look for candidate genes coding for enzymes belonging to the mevalonate biosynthetic pathway.
however, none of those genes were identified indicating a plausible partial loss of the mevalonate biosynthetic pathway.
as a consequence of this putative lost trait, mevalonate should be hypothetically imported from the environment , a situation that would suggest a new case of metabolic complementation between the bacterial endosymbiont and its host.
in favor of our hypothesis we point out that the mevalonate pathway plays a key role in insect metabolism as the precursor of juvenile hormone  and it is active in the fat body of b. germanica .
indirect evidence of the ability of mevalonate to diffuse and reach the endosymbiont is given by feeding experiments using mevalonate as a precursor of jh synthesis in the corpora allata .
moreover, the locus tag blbbge_ <dig> has been annotated as a homolog to ubie <cit> , which codes for an enzyme, a c-methyltransferase, that catalyzes reactions in both ubiquinone  and menaquinone  biosynthesis
<cit> .
in q biosynthesis, ubie catalyzes the conversion of 2-octaprenyl-6-methoxy- <dig> -benzoquinone to 2-octaprenyl-3-methyl-6-methoxy- <dig> -benzoquinone .
in mk biosynthesis, ubie catalyzes the conversion of demethylmenaquinone to menaquinone .
while the genome of blattabacterium has putative genes that code for the remaining activities for the mk biosynthesis, it has not any other annotated gene that accounts for the activities of the q biosynthesis.
as a consequence the activity ec  <dig> .
<dig>  has no biological meaning in the metabolic network of b. cuenoti and for that reason it should not be included in the model.
assuming a case of metabolic complementation between the bacteria and its host where mevalonate is a metabolite supplied by the host, a transport flux that allow the uptake of mevalonate by the cell was incorporated to the model.
after so um where recomputed just to find out that um <dig> was still non-functional.
a new inspection of the graph showed that there were two "dead end" reactions  for which no products were included as rnc metabolites.
a closer inspection of these reactions showed the following: reaction ec  <dig> .
<dig>  produces undecaprenyl diphosphate, a metabolite that works as a coenzyme in the biosynthesis of murein, whereas reaction ec  <dig> .
<dig>  is the last step in the biosynthesis of menaquinol, which is a known coenzyme that operates as an electron carrier.
as it is explained in section methods , both metabolites were then included in the biomass equation.
after so, um <dig> was completely solved, meaning that all reactions got unblocked .
um <dig> - nucleotide salvage pathway
in this case, the set of  <dig> reactions contained in the um belongs to the nucleotide salvage pathway.
the set of genes coding for these activities was analyzed and it was found that  <dig> of the reactions are orphans, i.e.
they do not have an associated coding gene .
the remaining  <dig> reactions were grouped according to its coding gene and this showed that  <dig> genes code these activities in the following way: seven reactions, that correspond to the activity  <dig> .
<dig> , are assigned to gene blbbge_ <dig> , annotated as purine-nucleoside phosphorylase; two reactions defined by the activities ec  <dig> .
<dig>  are associated to gene blbbge_ <dig> ; finally, the activity ec  <dig> .
<dig>  is assigned to the gene blbbge_ <dig> .
due to the great number of orphan reactions contained in um <dig> and based on the fact that all of these reactions are predicted as blocked, the first step to analyze this um was to remove its orphan reactions.
after so, the set of ums were recalculated to evaluate the impact of these changes.
it was found that the um splits into two ums: one of them formed by all the reactions associated with the activities ec  <dig> .
<dig>  together with the two reactions associated with activity ec  <dig> .
<dig> ; the other one was composed by the isolated reaction ec  <dig> .
<dig> .
in both situations, most of the metabolites involved in the ums were classified as rnp or rnc .
<dig>  however, in this case the reactions with no gene association , or wrong ec assignations, are represented with rounded rectangles and highlighted in yellow.
with the purpose of evaluating the functional assignment of these three genes, a close inspection of the genome annotation was done.
first, it was found that the annotation of blbbge_ <dig>  does not have strong evidence supporting that this gene code for the activity ec  <dig> .
<dig> .
thus, considering the lack of information supporting the association between the gene and the activity and taking into account that the model predicts that both reactions associated with activity ec  <dig> .
<dig>  are blocked, it is more plausible to assume that the activity does not take place in the metabolism of the endosymbiont.
second, the association of the enzyme activities ec  <dig> .
<dig>  and ec  <dig> .
<dig>  to genes blbbge_ <dig> and blbbge_ <dig> , respectively, is supported by the genome annotation.
however, the activity ec  <dig> .
<dig>  appears to have a broad substrate specificity and seven different reactions associated with this activity are included in model icg <dig>  all of them predicted to be blocked because their direct connection to an rnp, an rnc or to both.
in the case of the activity ec  <dig> .
<dig>  the reaction becomes isolated being one of its substrates and one of its products a rnp and a rnc metabolites, respectively.
due to the complexity of this situation and in the absence of experimental evidence that could help to solve these metabolic puzzles it was not possible to find the role of these activities in the metabolism of the bacterium, and thus the set of genes with its associated activities were excluded from the model until new experimental data shed some light over this problem.
um <dig> - pyridoxal 5-phosphate biosynthetic pathway
the seven reactions involved in this sub-graph correspond to the complete pyridoxal-5-phosphate pathway .
pyridoxal-5-phosphate, also called vitamin b <dig>  functions as a cofactor of different enzymes involved, among others, in transamination reactions required for the synthesis and catabolism of amino acids.
due to its importance in metabolism, this metabolite was included in the biomass equation in the same way that it has been done with other cofactors and coenzymes.
as a result all the reactions in the um become unblocked.
it is worth to note that animals, in particular insects, do not possess any biosynthetic pathway for pyridoxal-5-phosphate, and for this reason they need to take it from their diet in order to survive
<cit> .
as a consequence b. cuenoti may provide its host with this metabolite, suggesting another case of possible metabolic complementation.
um <dig> - lipopolysaccharide biosynthetic pathway
the um <dig> is composed by a linear chain of four reactions that are related to the biosynthetic pathways of different membrane lipids.
the first two reactions are involved in the palmitate biosynthetic pathway.
however, these reactions are not biochemically defined, but they are in turn the condensation of a set of reactions.
for example, the first reaction labeled as c120sn is the net sum of  <dig> activities that produce dodecanoyl-acp from acetoacetyl-acp.
moreover, the reaction kas <dig> is also the condensation of the activities ec  <dig> .
<dig>  and ec  <dig> .
<dig>  .
the other two reactions of the um are the first and second steps of the lipid iva biosynthetic pathway.
um <dig> was reformulated by decomposing the reactions c120sn and kas <dig> in its corresponding activities.
after so, the structure of the um was analyzed to find the following.
first, the activities ec  <dig> .
<dig>  and ec  <dig> .
<dig> , seem to be orphan in icg <dig>  second, the activities assigned to the lipid iva biosynthetic pathway were all orphan except for the activity ec  <dig> .
<dig> , which was assigned to the gene blbbge_ <dig> .
this scenario suggests that this pathway may be absent, and thus it could be the consequences of an error in the annotation of blbbge_ <dig> .
indeed, it was found that the set of ortholog genes identified in the genome of other sequenced genomes from diverse b. cuenoti strains have been annotated as coding for the activities ec  <dig> .
<dig>  and ec  <dig> .
<dig> .
as a consequence, the annotated activity of blbbge_ <dig>  was changed from ec  <dig> .
<dig>  to ec  <dig> .
<dig>  and ec  <dig> .
<dig> .
the remaining orphan activity present in the lipid iva pathway was also removed from the model because it is assumed that this pathway is not present in the metabolism of b. cuenoti.
um <dig> - siroheme biosynthetic pathway
a linear chain of four reactions composes the um <dig>  where the last reaction produces uroporphyrinogen iii, which was found to be a rnc .
this metabolite acts as substrate in the biosynthesis of siroheme, a prosthetic group which catalyzes the reduction of sulfite to sulfide and of nitrite to ammonia in the assimilation and dissimilation of sulfur and nitrogen compounds
<cit> .
after an inspection of b. cuenoti genome annotation, it was found that all the coding genes of the biosynthetic pathway of siroheme are present.
however, they had not been previously included in the metabolic model.
the pathway consists in four reactions  coded by two genes: blbbge_ <dig>  that codes the enzyme with activity ec  <dig> .
<dig>  and the gene blbbge_ <dig>  which codes an enzyme that catalyzes activities ec  <dig> .
<dig>  and ec  <dig> .
<dig> .
since siroheme is an important cofactor involved in sulfur and nitrogen metabolism, the cell would have to be able to maintain certain pool of this cofactor.
hence, during the bacterial growth phase, the organism will need some production of siroheme.
in order to take into account this fact, siroheme was included into the biomass equation.
after including this modification to the model, fba was applied to find a metabolic state that maximizes biomass production.
as expected, it was found that the four reactions included in um <dig> and the four new reactions added to the model showed non-zero flux under optimal state .
um <dig> - proline biosynthetic pathway
in this case two reactions were found: n-acetylornithine deacetylase coded by gene blbbg_ <dig>  and l-glutamate 5-semialdehyde dehydratase annotated to be spontaneous.
both reactions were annotated as belonging to the proline biosynthetic pathway.
in this scenario suggest that the pathway only lacks the last reaction step in order to be able to produce proline .
however a closer look to the annotation of gene blbbg_ <dig> shows that it has been assigned to code the following  <dig> activities:
<dig>  n-acetyl-l-glutamate 5-semialdehyde + h2o → acetate + l-glutamate 5-semialdehyde
<dig>  n-acetyl-l-ornithine + h2o → l-ornithine + acetate
<dig>  n-succinyl-l,l- <dig> -diaminopimelate + h2o → l,l-diaminopimelate + succinate
reaction  could not be found in the brenda database.
however it was found as an entry in bigg database
<cit> .
for these reason it is not clear whereas the reaction has been biochemically characterized or not.
in addition, experimental results show that proline is the most abundant amino acid in the cockroach’s hemolymph
<cit>  supporting the hypothesis that proline is provided by the host.
as a consequence, considering um <dig> as an annotation error seems a more plausible hypothesis and thus these activities were not included in the new version of the metabolic model.
um7–um <dig> - the case of isolated reactions
in the cases of um <dig> and um <dig>  both of them correspond to transport reactions associated to two different ions: fe2+ and mg2+ and, the associated exchange fluxes.
both ions are included in the biomass equation as described by other authors
<cit>  and in this way both ums got unblocked .
the isolated reaction that defines um <dig> is superoxide dismutase , a reaction that together with the activity ec  <dig> .
<dig>  works as a detoxification pathway against free radicals such as superoxide and hydrogen peroxide .
both activities have an associated coding gene and there is also experimental information suggesting the aerobic character of b. cuenoti <cit> .
taking together these facts it is expected that superoxide dismutase plays an important role in the metabolism of the endosymbiont.
the graph analysis showed the superoxide as a rnp metabolite.
the explanation of the previous finding rely in the fact that the processes of free radical formation, e.g.
as by product of aerobic respiration, is out of the scope of the model and hence the model does not include any reaction producing superoxide.
um <dig> is also a case of an isolated reaction  which catalyzes the activation of the apoprotein into acyl-carrier-protein , been this product a highly conserved carrier of acyl intermediates, important for fatty acid synthesis .
the process of protein biosynthesis is not included in the metabolic model of b. cuenoti  and for that reason the model doesn’t include any reaction producing the apoprotein.
as a consequence, the apoprotein becomes a rnp metabolite that blocks the activity ec  <dig> .
<dig> .
as in the case of um <dig>  some metabolite  is produced by a metabolic process that is out of the scope of the model and thus appears as a dead-end.
model update
the curation process described in the previous section resulted in the removal of  <dig> genes associated to  <dig> reactions from icg <dig> and the addition of  <dig> new genes corresponding to a total of  <dig> reactions.
thus, the new model version has two more genes  and for this reason has been named as imp <dig>  moreover, the reassignment and the inclusion of activities as well as the removal of orphan activities lead to elimination of  <dig> reactions  from icg <dig> and the inclusion of  <dig> new reactions  in imp <dig> .
since no experimental data was available for blattabacterium bge, chemical composition of e. coli, adapted from
<cit> , was used.
in particular, the stoichiometric coefficients of the new cofactor included to the biomass of imp <dig>  were those found in e. coli model ijo <dig>  even if blattabacterium and e. coli are phylogenetically very distant organisms, it is worth to note that these coefficients are approximations of the order of magnitude meant to capture the needs of an organism during growth in a qualitatively fashion.
for a detailed comparison between both models see additional file
<dig>
comparative analysis of the reaction subsets
the curation process involved addition and removal of reaction and metabolites as well as changes in the formulation of the biomass equation.
these changes affected the stoichiometric matrix, and then resulted in different structural properties of the network.
in particular, we have analyzed the organization of the reaction subsets .
table
<dig> summarizes the number of rs identified for the two models, as well the number of reaction within each rs.
the jaccard index was used as a measure of similarity between the rs from the two models.
this index was calculated for each pair of rs as the cardinality of the intersection over the cardinality of the union between both rss.
thus its value is bounded between  <dig> and  <dig>  the higher the index value, the higher is the number of reaction shared by both rs.
imp240
icg238
<dig>
.75
the first column indicates the number of reactions that belong to a rs.
the second and third columns indicate how many rs with a given number of reaction are in model imp <dig> and icg <dig>  respectively.
the fourth column shows how many rs pairs are equal in both models and for a given number of reactions.
in the last columns, the number out of the parenthesis indicate how many rs with a jaccard index greater than  <dig>  have been found, taking as the reference point the model imp <dig>  the number in parenthesis indicates the difference of reactions between a pair of rs in the following way: a positive number means that the rs in icg <dig> has this number of additional reactions whereas a negative indicate the opposite.
the major difference found in the reorganization was the presence of a rs composed by  <dig> reactions present in imp <dig> but not found in icg <dig>  this difference is due to the fact that model icg <dig> has all these  <dig> reactions condensed in a single step .
this is also the case of both the rs that contains  <dig> reactions and one of the rs composed by  <dig> reactions founded in imp <dig>  moreover, there is another rs of  <dig> reactions only present in imp <dig> that corresponds to the pyridoxal biosynthetic pathway.
this pathway was blocked in icg <dig>  and hence cannot be detected as an rs.
additionally, there are three rs almost equal in both models, but differing in only one reaction.
for example, the rs of  <dig> reactions in imp <dig> corresponds to the rs present in icg <dig> that contains  <dig> reactions.
despite the differences described above, no major changes were found in the organizations of rs between both models.
comparison of minimal medium
the in-silico minimal set of compounds needed for the endosymbiont in order to produce all biomass components was predicted for both model versions , using an optimization algorithm .
table
<dig> shows the results.
as it can be seen, all the metabolites included in the minimal medium predicted for model icg <dig> are included in the minimal medium predicted for imp <dig>  except for -dihydroorotate.
this metabolite is a precursor in the biosynthesis of pyrimidines.
the reason for this difference is that icg <dig> did not include the activity ec  <dig> .
<dig> , which catalyzes the conversion of n-carbamoyl-l-aspartate into -dihydroorotate, and thus the model predicts that -dihydroorotate should be uptaken from an external source.
the activity ec  <dig> .
<dig>  has been found to be coded by the gene blbbge_ <dig> , and its inclusion in the new model predicts that -dihydroorotate can be produced by the metabolism of the endosymbiont.
imp240
icg238
comparative table showing the in-silico predicted minimal medium for the new version model imp <dig> and previous version icg <dig>
the predictions using the model imp <dig> also suggest six new compounds that need to be present in the medium in order for the organism to be able to grow.
these sets of metabolites include: i) mevalonate, which is needed to synthesize menaquinol and 2-demetyl-menaquinol; ii) glycerol, which is phosphorylated by the activity ec  <dig> .
<dig> , and used in the biosynthesis of phosphatidylglycerol species.
in model igc <dig> glycerol-3-phosphate is produced from dihydroxyacetonephosphate through the reaction ec  <dig> .
<dig>  operating in reverse sense.
however, it was not possible to found any evidence suggesting that the reaction ec  <dig> .
<dig>  could operate in a reversible manner.
hence, if the reaction is considered as irreversible, then the cell cannot produce glycerol-3-phosphate.
then glycerol should be uptaken from the medium and phosphorylated inside the cell.
iii) porphobilinogen that is converted into hydroxymethylbilane, the first precursor in the synthesis of siroheme, which in turn requires fe2+; iv) the case of k+ and mg2+ are trivial: these ions become required in the minimal medium after their inclusion into the biomass equation.
network fragility
the set of essential metabolic genes were computed by an in-silico analysis of the model imp <dig>  these results were then compared with the set of essential genes previously predicted for model icg <dig> <cit> .
figure
<dig> summarizes the results of the comparative analysis between both models.
the set of  <dig> genes  predicted as essential in icg <dig> remains as essential in imp <dig>  in addition, a set of  <dig> genes was predicted as essential in imp <dig>  resulting in total of  <dig> metabolic genes predicted as essential .
as it is depicted in figure
<dig>   <dig> of those  <dig> genes were present in the previous version of the model , while the remaining  <dig> are new genes.
this increase in the predicted network fragility obtained after model curation is explained, in greater extent, due to the addition of new components to the biomass equation.
for example, the genes coding for the biosynthetic pathway of menaquinol where predicted as non-essential in model icg <dig> because the cofactor was not included in the biomass equation.
however, the cofactor must be essential to the organism, due to its strict aerobiosis, and then it should be included as a biomass component.
as a consequence most of the genes coding the biosynthetic pathway of menaquinol are predicted as essential in the new scenario.
moreover, the  <dig> new genes included in imp <dig> that are predicted as essentials, include two genes that code for three steps in the biosynthesis of siroheme, two genes that code for the transport of nicotinamide and mg+ <dig>  and the gene that encode ftsl, an essential cell division protein.
CONCLUSIONS
in this paper we have introduced a general definition of gap metabolite that allows its detection in a straightforward manner, even for the cases of upstream-non-produced and downstream-non-consumed metabolites.
moreover, a method for the detection of unconnected modules , defined as isolated sets of blocked reactions connected through gap metabolites have been proposed.
the visual representation of um can shed useful information that may help the model’s curator to solve inconsistencies.
furthermore, the present approach can be combined with existing tools in order to find a set of model modifications that solves the inconsistencies and thus improves model’s predictions.
this method was applied to the curation of the gsm of the cockroach endosymbiont b. cuenoti icg <dig>  in this way every blocked reaction detected in the model was successfully unblocked or alternatively removed, in those cases where there was not enough information supporting the existence of such reactions.
as an example, those reactions found as blocked and with no gene association were excluded.
moreover, new reactions were added to the model based on the careful revision of the genome annotation that allows the identification of gene functions previously not included, as well as the incorporation of new compounds into the biomass equation.
as a consequence of model curation a new gsm version of b. cuenoti, named as imp <dig>  is proposed.
the impact of these modifications, with respect to some structural properties of the networks, was analyzed by performing different in-silico analysis over each model’s version.
as a final commentary, the method here presented can be considered as a semi-automatic approach that has the advantage of allowing a quick representation of the gaps of the model but that needs the supervision of an expert in the biology of the studied organism.
this issue may be seen as a drawback of the method but, as in the case of automatic genome annotation, there is a trade-off between the degree of automation of the metabolic reconstruction and the quality of the generated model.
abbreviations
rnp: root-non-produced; rnc: root-non-consumed; dnp: downstream-non-produced; unc: upstream-non-consume; um: unconnected module; fba: flux balance analysis; cbm: constraint based modeling; gsm: genome scale model; rs: reactions subset; lp: linear programming; milp: mixed integer linear programming; gpr: gene-protein-reaction; ec: enzyme commission number; tc: transport commission number.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
jp, fm and mpl presented the original idea of the work and designed the study.
mpl conceived the algorithmic procedure, implemented the code, and carried out the in-silico experiments.
all authors contributed in designing research, analyzing the data, and writing the paper.
all authors have read and approved the manuscript.
supplementary material
additional file 1
schematic representation of ums.
the file contains schemes and a brief description of each um.
click here for file
additional file 2
new gsm model of b. cuenoti imp <dig>  spreadsheet with the model description.
the file contains three sheets corresponding to metabolites, reactions, and exchange fluxes, respectively.
click here for file
additional file 3
comparative table icg <dig> vs imp <dig>  spreadsheet describing the main differences between both model versions.
the file contains three sheets: the first one includes the genes added as well the genes that were removed; the second one presents the removed reactions; and the third one shows the set of added reactions.
click here for file
acknowledgments
mpl would like to thank fundación parque científico de madrid for allowing him to carry out the research activities corresponding to this paper.
financial support from spanish government is grateful acknowledged .

 metabolic reconstruction is the computational-based process that aims to elucidate the network of metabolites interconnected through reactions catalyzed by activities assigned to one or more genes. the reconstruction process begins with the functionally annotated genome of an organism. then, the identification of those genes whose putative products catalyze some biochemical reaction, i.e. gene products with an assigned enzyme commission number or transport commission number, should be done. this relational information can be organized in the so-called gene-protein-reaction association tables. in a further step gpr tables will be used to infer candidate metabolic pathways coded in the organism’s genome. in order to automate the process of reconstruction of a target organism’s metabolic network, computational methods have been previously developed that will yield a first draft. this first draft can be used to formulate a mathematical representation of an organism’s metabolism, termed as genome-scale model. the constraint-based-modeling is an approach that combines the stoichiometric analysis with optimization techniques to study genome-scale models. the cbm has been successfully used to predict metabolic capabilities such as growth rates, as well as systems responses to environmental or genetic perturbations. when applying cbm to an initial draft of a metabolic model, it is usual to find inconsistencies that can have different causes. in the initial stages of a metabolic reconstruction, due to annotation errors, as well as the existences of unknown enzyme functionality, gpr associations can be incorrectly established. thus, some reactions may be not included in the model draft. as a consequence, some metabolic pathways will contain gaps that will create dead-end metabolites. these metabolites appear in the model as only produced or only consumed by reactions, and hence will never reach a steady state different than the trivial, and then they will never participate in a feasible solution. they will in turn block any reaction in which they are involved. there are two classes of dead-end metabolites: i) root-non-produced metabolites i.e. metabolites that are only consumed by system’s reactions, and ii) root-non-consumed that includes those metabolites that are only produced by the network but never consumed. detection of rnp and rnc can be conducted by simply scanning the rows of the stoichiometric matrix. however, the absence of flow through metabolites rnp could be propagated downstream by blocking reactions and thus, additional metabolites would become gaps. those metabolites that become a gap as a consequence of some rnp metabolite are termed downstream-non-produced. in a symmetric way, upstream-non-consumed metabolites are defined as those metabolites that became a gap as a consequence of some present rnc metabolite. in general, the detection of dead-end metabolites and blocked reactions is referred commonly as the gap finding problem. model’s gaps can be solved by adding one or more reactions that allow connecting a dead-end metabolite with other metabolites of the network, a process known as gap-filling. in some cases, the incorporated reactions can be mapped into some coding gene. however, there are some situations where even if a set of candidate reactions that fill the gaps has been successfully predicted, it could not be possible to find the genes that code for these activities. in such cases the reactions are called orphan reactions. methods to predict candidate genes to be assigned to orphan reactions have also been previously developed. thus, the reconstruction of metabolic models is an iterative process in which the cbm plays an important role to detect inconsistencies that should be resolved or curated in order to improve model formulation. automated methods for model curation have been previously developed. in order to solve the gap-filling problem, an optimization-based method to identify the minimum number of reaction to be included in the model has been proposed by different authors. these methods rely in the use of mixed integer linear programming combined with universal reaction databases such as kegg, bigg or metacyc. other proposed approaches are based on the use of experimental information to detect inconsistencies with model predictions that may suggest errors in model formulation. even though automated methods for metabolic network curation are of an undoubted help to improve model formulation, there may be situations where the manual inspections of a curator are still needed. this is certainly the case of the reconstruction of networks from genomes that suffer reductive evolution and code for minimized metabolisms. during the establishment of symbiosis, metabolic redundancies with the host can result in the loss of enzymatic steps in the endosymbiont network, leading to the emergence of obligate metabolic complementation. these shared metabolic abilities take the form of interrupted pathways when the endosymbiotic network is reconstructed. thus, the problem of gap-filling is a complex decision making process where a visual representation of the inconsistencies can help to a model’s curator to understand how gap metabolites and blocked reactions are related and thus find the nature of these inconsistencies. in this paper we present a method that combines the cbm with an algorithm to compute connected components over bipartite-graphs. the presented method allows the detection of the isolated sets of blocked and gap metabolites and the way in which these are interconnected in what we have termed unconnected modules. then, the analysis of each individual unconnected module simplifies and clarifies the visual representation, and hence can be used to decide how gaps should be filled during the curation of the model. the availability of accurate gsm is especially relevant in the case of bacterial obligate endosymbionts since it is not possible the culture of these microorganisms and hence there is an intrinsic difficulty for obtaining direct experimental data from the system. in this case, modeling can serve as an appropriate proxy for functional characterizations. for instance, metabolic modeling has been successfully used in the case of buchnera aphidicola to evaluate the role of the endosymbiont and the host in nitrogen metabolism, in sodalis glossinidius to characterize intermediate steps during the reduction of the network as a result of the interaction with the host metabolism, and in blattabacterium cuenoti to better understand the striking conservation of the endosymbiont metabolism along the evolutionary time as well as its putative role in the nitrogen economy of the system. in this paper we have used the gsm from b. cuenoti to test the proposed method of curation. our study allows upgrading icg to a more accurate model version named imp and it can be used as a guide for systems biology experimental explorations of the interaction between cockroaches and their endosymbionts. methods constraint-based modeling the study of the structural properties of biochemical reactions network relies on the analysis of the stoichiometric matrix. let denote by n the stoichiometric matrix associated to a certain metabolic network with m rows and n columns corresponding to the number of metabolites and reactions, respectively. in the following i and j will refer to the set of metabolite indexes and reaction indexes, respectively. moreover, the set of reaction indexes j will be partitioned into two disjoint subsets: the set jint which contains the indexes of internal fluxes, i.e. the biochemical reactions that take place inside the cell, as well the transport reactions that operate between the cell and the surrounding medium. on the other hand, the set jex contains the indexes of the exchange fluxes, which are the auxiliary variables used to represent the rate at which certain metabolites are consumed or produced by the system. there will be only one exchange flux per metabolite and these fluxes will be associated to the metabolites belonging to the extra-cellular compartment. by convention, the activity of the exchange fluxes is defined as positive or negative if the metabolite is produced or consumed by the system, respectively. the cbm approach relays on the use of different kinds of constraints represented by mathematical equations to define the so called flux space f, i.e. the set of all flux distributions compatible with the given constraints. the steady state condition is imposed over the mass balance equation of each metabolite of the network yielding the following homogeneous system of linear equations: n.v= where the vector v is a flux distribution compatible with the steady-state condition. moreover, lower and upper bounds are imposed over each reaction to represent additional constraints. for instance, the thermodynamic constraints that make some reactions to be irreversible are represented by setting to zero the lower bound. besides, the surrounding environment of a metabolic system can be modeled by setting bounds over the exchange fluxes. for example, if a given metabolite is available in the medium and thus can be consumed by the system, the lower bound of the corresponding exchange flux should have a negative value. the lower and upper bounds imposed over each flux can be written as the following systems of linear inequalities: vjlb≤vj≤vjub∀jεj where vj is the activity through the flux j, whereas vjlb and vjub are its lower and upper bounds, respectively. together, the homogeneous system of linear equations and the system of linear inequalities yields the mathematical representation of the flux space f, expressed as: f=vεrn:n.v= vjlb≤vj≤vjub∀jεj blocked reactions a reaction in a metabolic model is defined as blocked under a given medium condition if it cannot display a steady-state flux other than zero: jεjblocked⇔vj= ∀vεf where jblocked is the set of blocked reaction indexes. this set can be computed solving a set of linear programs, as proposed by burgard et al.. the approach consists in calculating the minimum and maximum flux value through each reaction of the system. when the maximum and minimum values found for a given reaction are both equal to zero, the reaction is said to be blocked under the defined medium condition. the formulation of the set of linear programs is the following: min/max:vj∀jεjs.t.n⋅v=0vjlb≤vj≤vjub∀jεj gap metabolites gap metabolites in gsm are defined as those vertexes of the network through which there can be no steady state flow. while the detection of rnp and rnc metabolites is straightforward by scanning of each row of the stoichiometric matrix n, the case of detecting unp and dnc metabolites cannot be accomplished by a simple inspection of the entries of n. however, based on the definition of blocked reaction we found a way to define gap metabolites that allow its identification in a straightforward manner. definition: a metabolite ί ∈ i in a network under steady-state is a gap if and only if all the reactions in which its participate belongs to the set of blocked reactions j blocked. this statement implies that there cannot exists a stationary flow through this metabolite. thus, if we name the set of reactions in which a metabolite ί participate as: σi=jεj:nij≠0∀iεi then, the set of gap metabolites igap ⊆ i can be defined as follows: i∈igap⇔σi⊆jblocked hence, the detection of gap metabolites can be accomplished by finding the set jblocked and applying equations and for each metabolite. the given definition for a gap metabolite doesn’t make any distinction between the different classes of gap defined. although there is no general procedure for the classification of gap metabolites as rnp, rnc, unp, or dnc once they have been found, this is possible in simple cases. for example, if a gap metabolite ί is involved only in an irreversible reaction by which it is consumed, or if all the reactions in which this metabolite is involved are irreversible and in all it is consumed, then metabolite will be a rnp. a similar reasoning can we argued for the symmetric case, leading to the identification of a rnc. for more complex cases, a visual inspection of the underling "unconnected module" may help to the classification. the coenzyme pseudo-gap problem in general gap metabolites can be identified using its relation to the set of blocked reactions jblocked as it was explained in previous subsection. however, there could be some special metabolites that are not "gap" under the definition given by. nevertheless these metabolites may be the cause that certain reactions get blocked, as the example depicted in figure  shows. the metabolite d is not a gap but the mass balances equation for d* implies that v is equal to v and the mass balance equation for d implies v + v = v as a consequence of these relations v is restricted to zero, i.e. v is a blocked reaction. this effect propagates downstream to v v and v this is why we name d as a pseudo-gap metabolite. in order to unblock v which in term will unblock v v and v there must be added a sink for d or d*. these kinds of situations may take place when the biosynthetic pathway of a coenzyme is present in a model, but there are no fluxes draining or degrading the coenzyme produced by this pathway. these metabolites may be involved in conserved moieties, and in such cases they will be consumed and regenerated in a cyclic manner. as a consequence they could not be detected as gap metabolites because they will participate in at least two active reactions. however, if the biosynthetic pathway for a coenzyme is included in a metabolic model, the net production of a coenzyme will not occur under steady state unless some flux consumes it. hence, the reactions involved in the biosynthesis pathway may become blocked. a common approach to solve this problem is to include the coenzyme-like metabolite into the biomass equation or alternatively to introduce an exchange flux that can drain the metabolite out of the system. either of these two situations is equivalent to the addition of the above-mentioned sink. the way to detect the pseudo-gap metabolites may be summarized as follow. when an um overlap with the biosynthetic pathway of a certain cofactor, two different situations can be found: the cofactor is included in the um as a gap or it is not. if it is a gap, the cofactor must be included in the biomass equation in order to solve the um. if the cofactor is not included in the um, then it is quite probable that the cofactor may be involved in a conservation relation which connect the cofactor to active reactions and for this reason is not detected as a gap under the definition given by equation. however, even not detected as gap, the cofactor could be the underlying cause that the um become blocked, and for this reason was termed the cofactor a "pseudo-gap" as has been pointed out at the beginning of this section. the conservation relations can be detected by the analysis of the conserved moieties, which were calculated as previously described. finally, by adding the corresponding cofactor to the biomass equation the um can be solved. identification of unconnected modules when analyzing the set of blocked reactions and gap metabolites of a metabolic model, it is common to find relations between both sets due to the fact that blocked reactions may be connected with other blocked reactions through gap metabolite. in some cases, blocked reactions are directly connected to a rnp/rnc metabolite. however, this may be not the case of other blocked reactions that may be not linked directly to a rnp or rnc metabolite. for example, as it can be seen in figure  1a, reaction v although blocked, is not directly connected to a rnp metabolite, but indirectly through metabolite through a set dnc metabolites and other blocked reactions. similarly occurs in the symmetric case. as a consequence of these relations, it is possible to systematically establish how the blocked reactions are connected through the gap metabolites. any metabolic network can be represented as a directed bipartite graph by considering the sets of vertex v = iuj. then, a directed edge or arc will exist between a metabolite and a reaction if the metabolite participates in the reaction. the direction of the arc will be incident to the reaction if the metabolite is a reactant and incident to the metabolite if it is a product. in the following, the graph associated to a metabolic network will be referred as the metabolic graph. once the metabolic graph is constructed it is possible to consider any possible sub-graph by selecting a pair of subsets i' ⊆ i and j' ⊆ j of metabolites and reactions respectively. in particular, it is possible to consider the sub-graph defined by the subset of vertex v = igapujblocked. this sub-graph will contain the relation that exists between gap metabolites and blocked reactions. moreover, the set of connected components can be computed over this graph. in this context, each connected component can be interpreted as a "module" of the metabolism that becomes inactive or unconnected, possibly as a consequence of model inconsistencies, such as the presences of a set of rnp/rnc metabolites. for this reason the set of connected component will be referred as an unconnected modules. in simple cases, the reason that causes a certain um to be unconnected from the rest of the network may be found by visual inspection of the graph representing the um. in such cases connectivity restoring of the rnc and rnp metabolites present in the um may, in general, solve the problem of the unc and dnp metabolites. the set of elementary operations that can be applied to solve ums has been discussed by kumar and collaborators, and it includes: addition of biochemical reactions of transport, incorporation of exchange fluxes and relaxation of irreversibility constraint of some reactions. more complex situations may include cases such as the pseudo-gap problem described in the previous section. flux balance analysis the flux balance analysis is an approach that combines the description of the flux space f defined by equation with optimization techniques to find a flux distribution v that maximizes the growth rate. this problem can formulated as a linear program and can be solved with standard techniques of linear programming. max:vbiomasss.t.n⋅v=0vjlb≤vj≤vjub∀jεj where the flux vbiomass represents the growth rate of the organism. this flux, also refereed as the biomass equation, includes all the metabolites that are biomass components, in its specific proportions. in-silico knockout experiments the fragility analysis of a network was performed by simulating knockout experiments for each metabolic gene included in a gsm. an in-silico knockout experiment for a given gene consists in bound to zero the flux for each reaction coded by the gene, which are inferred through the grp association table. after so, fba is used to find the maximal value of the biomass reaction under the genetic perturbation. if the optimal value of biomass reaction is lower than a certain threshold, the knockout is said to be lethal, otherwise the gene is not essential. the procedure is performed over all the genes in the model. minimal medium prediction the minimal medium may be defined as the smaller set of metabolites that should be present in the medium condition in order a feasible flux distribution v to exist, with a biomass production rate greater than zero. the minimal medium was calculated by solving a milp algorithm as previously described. the milp algorithm is the following: min:∑jεjexyjs.t.n⋅v=0vjlb≤vj≤vjub∀jεjintvbiomass≥vbiomasslbvjlbyj≤vj∀jεjexyjε ∀jεjex the algorithm requires the incorporation of a set of binary yj variables, one for each exchange flux j ε jex. moreover, a set of constraints that relates each binary variable with its corresponding exchanges flux should be incorporated to the problem. then, whenever a binary variable takes the zero value, the corresponding exchange flux is constrained to take a value greater or equal than zero. additionally, the lower bound corresponding to the biomass production flux should be set to a cutoff value greater than zero in order to guarantee a positive growth rate. finally, the optimization target is defined in such a way that minimizes the number of active exchange fluxes with a negative value. due to the fact that each exchange flux is related with one extracellular metabolite this is equivalent to find the minimal set of metabolites that the system must consume in order to produce a biomass flux greater than zero. detection of reaction subsets a reaction subsets or full coupling sets in a stoichiometric network is a group of reactions that operate together in fixed flux proportions for any flux distribution. due to the fact that the relations between enzymes and reactions are not always biunivocal, the reaction subset does not always match the concept of enzyme subset previously introduce by pfeiffer et al., and for this reason the term reaction subset seams more appropriate. the rss are structural invariant of network and for that they are independent on the kinetic parameters of the system. moreover, they shed valuable information that may help to understand how the network is regulated.for these reasons the concept rs is important for the analysis of metabolic networks. in order to compute reaction subset the following pre-processing steps were applied to the network. first, rows and columns corresponding to gap metabolites and blocked reaction respectively were removed from the stoichiometric matrix. as it was previously described, the constant biomass composition imposed by the stoichiometry of the biomass equation was relaxed by removing the corresponding column, while allowing each biomass component to be drained from the system in an independent way. then, the identification of rs was done by using the algorithm described in. computational tools constraint-based analysis was performed using the python-based toolbox cobrapy. lp and milp problems were solved using the gurobi solver acceded through cobrapy. identification of conservation relations was done by computation the set of extreme rays using the polco package. identification of rs was done by using an implementation based on python of the algorithm described in. the detection of the connected components of a graph was done using the implementation of the algorithm available in the igraph library. graphs were drawn using the yed graph editor. all computation was done on a desktop computer with an intel® core™ i cpu processor, with gib, running under fedora linux os.  the first step in the analysis of the gsm of b. cuenoti icg was to find the sets jblocked and igap of blocked reactions and gap metabolites, respectively. the results showed that reactions over a total of are blocked under any medium condition. using this information a set of metabolites over a total of were detected as gaps. a bipartite graph representation of the metabolic network was constructed, and the sub-graph defined by the subset of vertex igap u jblocked was selected. computation of connected components over this sub-graph allows identifying um. after so, each gap metabolite was classified in one of the four categories: rnp, rnc, unc and dnp. a description of each different um found in icg is summarized in table  which was sorted according to the number of reactions included in each um. um identified in b. cuenoti icg gsm. metabolite abbreviations: mev, 2ombzl, 2ommbl, hxan, xan, r1p, 2dr1p, ura 1pyr5c, 4hthr, u3hga -d-glucosamine), uppg, apoacp. the relation between an um and a subsystem was established according to the most frequent subsystem associated to the reactions of the um. when analyzing the reactions participation of each um, it was found that in many of the cases the set of reactions belonging to an um overlapped with known biochemical subsystems. for example, all the reaction included in um belongs to the menaquinol biosynthetic pathway. it is worth to note that in some cases ums could be composed by an isolated reaction as the cases of um and um after the identification of all um applying the proposed approach, each sub-graph was drawn independently. visual inspection of each graph was done to detect possible inconsistencies in the network formulation that may help to carry out the manual curation of a gsm. um - menaquinol biosynthetic pathway the biggest um found in the metabolic model was first analyzed. figure  shows the graph that represents the um. a visual inspection of these figure shows that um includes all the reactions and metabolite that conforms the menaquinol biosynthetic pathway. the figure also shows that from the gap metabolites included in the um, two are rnp metabolites, and one is an rnc metabolite, while the remaining metabolites are then upstream/downstream gap metabolites. mevalonate is a precursor for the biosynthesis of isopentenyl diphosphate that in turn is involved in the biosynthetic pathway of many cofactors such as menaquinol and 2-demetyl-menaquinol. an inspection of the current genome annotation of b. cuenoti was done to look for candidate genes coding for enzymes belonging to the mevalonate biosynthetic pathway. however, none of those genes were identified indicating a plausible partial loss of the mevalonate biosynthetic pathway. as a consequence of this putative lost trait, mevalonate should be hypothetically imported from the environment, a situation that would suggest a new case of metabolic complementation between the bacterial endosymbiont and its host. in favor of our hypothesis we point out that the mevalonate pathway plays a key role in insect metabolism as the precursor of juvenile hormone and it is active in the fat body of b. germanica. indirect evidence of the ability of mevalonate to diffuse and reach the endosymbiont is given by feeding experiments using mevalonate as a precursor of jh synthesis in the corpora allata. moreover, the locus tag blbbge_ has been annotated as a homolog to ubie, which codes for an enzyme, a c-methyltransferase, that catalyzes reactions in both ubiquinone and menaquinone biosynthesis. in q biosynthesis, ubie catalyzes the conversion of 2-octaprenyl-6-methoxy- -benzoquinone to 2-octaprenyl-3-methyl-6-methoxy- -benzoquinone. in mk biosynthesis, ubie catalyzes the conversion of demethylmenaquinone to menaquinone. while the genome of blattabacterium has putative genes that code for the remaining activities for the mk biosynthesis, it has not any other annotated gene that accounts for the activities of the q biosynthesis. as a consequence the activity ec. has no biological meaning in the metabolic network of b. cuenoti and for that reason it should not be included in the model. assuming a case of metabolic complementation between the bacteria and its host where mevalonate is a metabolite supplied by the host, a transport flux that allow the uptake of mevalonate by the cell was incorporated to the model. after so um where recomputed just to find out that um was still non-functional. a new inspection of the graph showed that there were two "dead end" reactions for which no products were included as rnc metabolites. a closer inspection of these reactions showed the following: reaction ec. produces undecaprenyl diphosphate, a metabolite that works as a coenzyme in the biosynthesis of murein, whereas reaction ec. is the last step in the biosynthesis of menaquinol, which is a known coenzyme that operates as an electron carrier. as it is explained in section methods, both metabolites were then included in the biomass equation. after so, um was completely solved, meaning that all reactions got unblocked. um - nucleotide salvage pathway in this case, the set of reactions contained in the um belongs to the nucleotide salvage pathway. the set of genes coding for these activities was analyzed and it was found that of the reactions are orphans, i.e. they do not have an associated coding gene. the remaining reactions were grouped according to its coding gene and this showed that genes code these activities in the following way: seven reactions, that correspond to the activity., are assigned to gene blbbge_, annotated as purine-nucleoside phosphorylase; two reactions defined by the activities ec. are associated to gene blbbge_; finally, the activity ec. is assigned to the gene blbbge_. due to the great number of orphan reactions contained in um and based on the fact that all of these reactions are predicted as blocked, the first step to analyze this um was to remove its orphan reactions. after so, the set of ums were recalculated to evaluate the impact of these changes. it was found that the um splits into two ums: one of them formed by all the reactions associated with the activities ec. together with the two reactions associated with activity ec.; the other one was composed by the isolated reaction ec.. in both situations, most of the metabolites involved in the ums were classified as rnp or rnc. however, in this case the reactions with no gene association, or wrong ec assignations, are represented with rounded rectangles and highlighted in yellow. with the purpose of evaluating the functional assignment of these three genes, a close inspection of the genome annotation was done. first, it was found that the annotation of blbbge_ does not have strong evidence supporting that this gene code for the activity ec.. thus, considering the lack of information supporting the association between the gene and the activity and taking into account that the model predicts that both reactions associated with activity ec. are blocked, it is more plausible to assume that the activity does not take place in the metabolism of the endosymbiont. second, the association of the enzyme activities ec. and ec. to genes blbbge_ and blbbge_, respectively, is supported by the genome annotation. however, the activity ec. appears to have a broad substrate specificity and seven different reactions associated with this activity are included in model icg all of them predicted to be blocked because their direct connection to an rnp, an rnc or to both. in the case of the activity ec. the reaction becomes isolated being one of its substrates and one of its products a rnp and a rnc metabolites, respectively. due to the complexity of this situation and in the absence of experimental evidence that could help to solve these metabolic puzzles it was not possible to find the role of these activities in the metabolism of the bacterium, and thus the set of genes with its associated activities were excluded from the model until new experimental data shed some light over this problem. um - pyridoxal 5-phosphate biosynthetic pathway the seven reactions involved in this sub-graph correspond to the complete pyridoxal-5-phosphate pathway. pyridoxal-5-phosphate, also called vitamin b functions as a cofactor of different enzymes involved, among others, in transamination reactions required for the synthesis and catabolism of amino acids. due to its importance in metabolism, this metabolite was included in the biomass equation in the same way that it has been done with other cofactors and coenzymes. as a result all the reactions in the um become unblocked. it is worth to note that animals, in particular insects, do not possess any biosynthetic pathway for pyridoxal-5-phosphate, and for this reason they need to take it from their diet in order to survive. as a consequence b. cuenoti may provide its host with this metabolite, suggesting another case of possible metabolic complementation. um - lipopolysaccharide biosynthetic pathway the um is composed by a linear chain of four reactions that are related to the biosynthetic pathways of different membrane lipids. the first two reactions are involved in the palmitate biosynthetic pathway. however, these reactions are not biochemically defined, but they are in turn the condensation of a set of reactions. for example, the first reaction labeled as c120sn is the net sum of activities that produce dodecanoyl-acp from acetoacetyl-acp. moreover, the reaction kas is also the condensation of the activities ec. and ec.. the other two reactions of the um are the first and second steps of the lipid iva biosynthetic pathway. um was reformulated by decomposing the reactions c120sn and kas in its corresponding activities. after so, the structure of the um was analyzed to find the following. first, the activities ec. and ec., seem to be orphan in icg second, the activities assigned to the lipid iva biosynthetic pathway were all orphan except for the activity ec., which was assigned to the gene blbbge_. this scenario suggests that this pathway may be absent, and thus it could be the consequences of an error in the annotation of blbbge_. indeed, it was found that the set of ortholog genes identified in the genome of other sequenced genomes from diverse b. cuenoti strains have been annotated as coding for the activities ec. and ec.. as a consequence, the annotated activity of blbbge_ was changed from ec. to ec. and ec.. the remaining orphan activity present in the lipid iva pathway was also removed from the model because it is assumed that this pathway is not present in the metabolism of b. cuenoti. um - siroheme biosynthetic pathway a linear chain of four reactions composes the um where the last reaction produces uroporphyrinogen iii, which was found to be a rnc. this metabolite acts as substrate in the biosynthesis of siroheme, a prosthetic group which catalyzes the reduction of sulfite to sulfide and of nitrite to ammonia in the assimilation and dissimilation of sulfur and nitrogen compounds. after an inspection of b. cuenoti genome annotation, it was found that all the coding genes of the biosynthetic pathway of siroheme are present. however, they had not been previously included in the metabolic model. the pathway consists in four reactions coded by two genes: blbbge_ that codes the enzyme with activity ec. and the gene blbbge_ which codes an enzyme that catalyzes activities ec. and ec.. since siroheme is an important cofactor involved in sulfur and nitrogen metabolism, the cell would have to be able to maintain certain pool of this cofactor. hence, during the bacterial growth phase, the organism will need some production of siroheme. in order to take into account this fact, siroheme was included into the biomass equation. after including this modification to the model, fba was applied to find a metabolic state that maximizes biomass production. as expected, it was found that the four reactions included in um and the four new reactions added to the model showed non-zero flux under optimal state. um - proline biosynthetic pathway in this case two reactions were found: n-acetylornithine deacetylase coded by gene blbbg_ and l-glutamate 5-semialdehyde dehydratase annotated to be spontaneous. both reactions were annotated as belonging to the proline biosynthetic pathway. in this scenario suggest that the pathway only lacks the last reaction step in order to be able to produce proline. however a closer look to the annotation of gene blbbg_ shows that it has been assigned to code the following activities: n-acetyl-l-glutamate 5-semialdehyde + h2o → acetate + l-glutamate 5-semialdehyde n-acetyl-l-ornithine + h2o → l-ornithine + acetate n-succinyl-l,l- -diaminopimelate + h2o → l,l-diaminopimelate + succinate reaction could not be found in the brenda database. however it was found as an entry in bigg database. for these reason it is not clear whereas the reaction has been biochemically characterized or not. in addition, experimental results show that proline is the most abundant amino acid in the cockroach’s hemolymph supporting the hypothesis that proline is provided by the host. as a consequence, considering um as an annotation error seems a more plausible hypothesis and thus these activities were not included in the new version of the metabolic model. um7–um - the case of isolated reactions in the cases of um and um both of them correspond to transport reactions associated to two different ions: fe2+ and mg2+ and, the associated exchange fluxes. both ions are included in the biomass equation as described by other authors and in this way both ums got unblocked. the isolated reaction that defines um is superoxide dismutase, a reaction that together with the activity ec. works as a detoxification pathway against free radicals such as superoxide and hydrogen peroxide. both activities have an associated coding gene and there is also experimental information suggesting the aerobic character of b. cuenoti. taking together these facts it is expected that superoxide dismutase plays an important role in the metabolism of the endosymbiont. the graph analysis showed the superoxide as a rnp metabolite. the explanation of the previous finding rely in the fact that the processes of free radical formation, e.g. as by product of aerobic respiration, is out of the scope of the model and hence the model does not include any reaction producing superoxide. um is also a case of an isolated reaction which catalyzes the activation of the apoprotein into acyl-carrier-protein, been this product a highly conserved carrier of acyl intermediates, important for fatty acid synthesis. the process of protein biosynthesis is not included in the metabolic model of b. cuenoti and for that reason the model doesn’t include any reaction producing the apoprotein. as a consequence, the apoprotein becomes a rnp metabolite that blocks the activity ec.. as in the case of um some metabolite is produced by a metabolic process that is out of the scope of the model and thus appears as a dead-end. model update the curation process described in the previous section resulted in the removal of genes associated to reactions from icg and the addition of new genes corresponding to a total of reactions. thus, the new model version has two more genes and for this reason has been named as imp moreover, the reassignment and the inclusion of activities as well as the removal of orphan activities lead to elimination of reactions from icg and the inclusion of new reactions in imp. since no experimental data was available for blattabacterium bge, chemical composition of e. coli, adapted from, was used. in particular, the stoichiometric coefficients of the new cofactor included to the biomass of imp were those found in e. coli model ijo even if blattabacterium and e. coli are phylogenetically very distant organisms, it is worth to note that these coefficients are approximations of the order of magnitude meant to capture the needs of an organism during growth in a qualitatively fashion. for a detailed comparison between both models see additional file comparative analysis of the reaction subsets the curation process involved addition and removal of reaction and metabolites as well as changes in the formulation of the biomass equation. these changes affected the stoichiometric matrix, and then resulted in different structural properties of the network. in particular, we have analyzed the organization of the reaction subsets. table  summarizes the number of rs identified for the two models, as well the number of reaction within each rs. the jaccard index was used as a measure of similarity between the rs from the two models. this index was calculated for each pair of rs as the cardinality of the intersection over the cardinality of the union between both rss. thus its value is bounded between and the higher the index value, the higher is the number of reaction shared by both rs. imp240 icg238.75 the first column indicates the number of reactions that belong to a rs. the second and third columns indicate how many rs with a given number of reaction are in model imp and icg respectively. the fourth column shows how many rs pairs are equal in both models and for a given number of reactions. in the last columns, the number out of the parenthesis indicate how many rs with a jaccard index greater than have been found, taking as the reference point the model imp the number in parenthesis indicates the difference of reactions between a pair of rs in the following way: a positive number means that the rs in icg has this number of additional reactions whereas a negative indicate the opposite. the major difference found in the reorganization was the presence of a rs composed by reactions present in imp but not found in icg this difference is due to the fact that model icg has all these reactions condensed in a single step. this is also the case of both the rs that contains reactions and one of the rs composed by reactions founded in imp moreover, there is another rs of reactions only present in imp that corresponds to the pyridoxal biosynthetic pathway. this pathway was blocked in icg and hence cannot be detected as an rs. additionally, there are three rs almost equal in both models, but differing in only one reaction. for example, the rs of reactions in imp corresponds to the rs present in icg that contains reactions. despite the differences described above, no major changes were found in the organizations of rs between both models. comparison of minimal medium the in-silico minimal set of compounds needed for the endosymbiont in order to produce all biomass components was predicted for both model versions, using an optimization algorithm. table  shows the results. as it can be seen, all the metabolites included in the minimal medium predicted for model icg are included in the minimal medium predicted for imp except for -dihydroorotate. this metabolite is a precursor in the biosynthesis of pyrimidines. the reason for this difference is that icg did not include the activity ec., which catalyzes the conversion of n-carbamoyl-l-aspartate into -dihydroorotate, and thus the model predicts that -dihydroorotate should be uptaken from an external source. the activity ec. has been found to be coded by the gene blbbge_, and its inclusion in the new model predicts that -dihydroorotate can be produced by the metabolism of the endosymbiont. imp240 icg238 comparative table showing the in-silico predicted minimal medium for the new version model imp and previous version icg the predictions using the model imp also suggest six new compounds that need to be present in the medium in order for the organism to be able to grow. these sets of metabolites include: i) mevalonate, which is needed to synthesize menaquinol and 2-demetyl-menaquinol; ii) glycerol, which is phosphorylated by the activity ec., and used in the biosynthesis of phosphatidylglycerol species. in model igc glycerol-3-phosphate is produced from dihydroxyacetonephosphate through the reaction ec. operating in reverse sense. however, it was not possible to found any evidence suggesting that the reaction ec. could operate in a reversible manner. hence, if the reaction is considered as irreversible, then the cell cannot produce glycerol-3-phosphate. then glycerol should be uptaken from the medium and phosphorylated inside the cell. iii) porphobilinogen that is converted into hydroxymethylbilane, the first precursor in the synthesis of siroheme, which in turn requires fe2+; iv) the case of k+ and mg2+ are trivial: these ions become required in the minimal medium after their inclusion into the biomass equation. network fragility the set of essential metabolic genes were computed by an in-silico analysis of the model imp these results were then compared with the set of essential genes previously predicted for model icg. figure  summarizes the results of the comparative analysis between both models. the set of genes predicted as essential in icg remains as essential in imp in addition, a set of genes was predicted as essential in imp resulting in total of metabolic genes predicted as essential. as it is depicted in figure  of those genes were present in the previous version of the model, while the remaining are new genes. this increase in the predicted network fragility obtained after model curation is explained, in greater extent, due to the addition of new components to the biomass equation. for example, the genes coding for the biosynthetic pathway of menaquinol where predicted as non-essential in model icg because the cofactor was not included in the biomass equation. however, the cofactor must be essential to the organism, due to its strict aerobiosis, and then it should be included as a biomass component. as a consequence most of the genes coding the biosynthetic pathway of menaquinol are predicted as essential in the new scenario. moreover, the new genes included in imp that are predicted as essentials, include two genes that code for three steps in the biosynthesis of siroheme, two genes that code for the transport of nicotinamide and mg+ and the gene that encode ftsl, an essential cell division protein.  in this paper we have introduced a general definition of gap metabolite that allows its detection in a straightforward manner, even for the cases of upstream-non-produced and downstream-non-consumed metabolites. moreover, a method for the detection of unconnected modules, defined as isolated sets of blocked reactions connected through gap metabolites have been proposed. the visual representation of um can shed useful information that may help the model’s curator to solve inconsistencies. furthermore, the present approach can be combined with existing tools in order to find a set of model modifications that solves the inconsistencies and thus improves model’s predictions. this method was applied to the curation of the gsm of the cockroach endosymbiont b. cuenoti icg in this way every blocked reaction detected in the model was successfully unblocked or alternatively removed, in those cases where there was not enough information supporting the existence of such reactions. as an example, those reactions found as blocked and with no gene association were excluded. moreover, new reactions were added to the model based on the careful revision of the genome annotation that allows the identification of gene functions previously not included, as well as the incorporation of new compounds into the biomass equation. as a consequence of model curation a new gsm version of b. cuenoti, named as imp is proposed. the impact of these modifications, with respect to some structural properties of the networks, was analyzed by performing different in-silico analysis over each model’s version. as a final commentary, the method here presented can be considered as a semi-automatic approach that has the advantage of allowing a quick representation of the gaps of the model but that needs the supervision of an expert in the biology of the studied organism. this issue may be seen as a drawback of the method but, as in the case of automatic genome annotation, there is a trade-off between the degree of automation of the metabolic reconstruction and the quality of the generated model. abbreviations rnp: root-non-produced; rnc: root-non-consumed; dnp: downstream-non-produced; unc: upstream-non-consume; um: unconnected module; fba: flux balance analysis; cbm: constraint based modeling; gsm: genome scale model; rs: reactions subset; lp: linear programming; milp: mixed integer linear programming; gpr: gene-protein-reaction; ec: enzyme commission number; tc: transport commission number. competing interests the authors declare that they have no competing interests. authors’ contributions jp, fm and mpl presented the original idea of the work and designed the study. mpl conceived the algorithmic procedure, implemented the code, and carried out the in-silico experiments. all authors contributed in designing research, analyzing the data, and writing the paper. all authors have read and approved the manuscript. supplementary material additional file 1 schematic representation of ums. the file contains schemes and a brief description of each um. click here for file additional file 2 new gsm model of b. cuenoti imp spreadsheet with the model description. the file contains three sheets corresponding to metabolites, reactions, and exchange fluxes, respectively. click here for file additional file 3 comparative table icg vs imp spreadsheet describing the main differences between both model versions. the file contains three sheets: the first one includes the genes added as well the genes that were removed; the second one presents the removed reactions; and the third one shows the set of added reactions. click here for file acknowledgments mpl would like to thank fundación parque científico de madrid for allowing him to carry out the research activities corresponding to this paper. financial support from spanish government is grateful acknowledged.
reconstructed models may contain inconsistencies that appear as gap metabolites and blocked reactions. moreover, a method for the detection of unconnected modules, defined as isolated sets of blocked reactions connected through gap metabolites, is proposed. metabolic reconstruction is the computational-based process that aims to elucidate the network of metabolites interconnected through reactions catalyzed by activities assigned to one or more genes. we introduce a general definition of gap metabolite that allows its detection in a straightforward manner. the outcome of its application to the genome-scale model b. cuenoti icg is a more accurate model version named as b. cuenoti imp although automatic methods for solving this problem have been previously developed, there are many situations where manual curation is still needed.
 metabolic reconstruction is the computational-based process that aims to elucidate the network of metabolites interconnected through reactions catalyzed by activities assigned to one or more genes. reconstructed models may contain inconsistencies that appear as gap metabolites and blocked reactions. although automatic methods for solving this problem have been previously developed, there are many situations where manual curation is still needed.  we introduce a general definition of gap metabolite that allows its detection in a straightforward manner. moreover, a method for the detection of unconnected modules, defined as isolated sets of blocked reactions connected through gap metabolites, is proposed. the method has been successfully applied to the curation of icg the genome-scale metabolic model for the bacterium blattabacterium cuenoti, obligate endosymbiont of cockroaches.  we found the proposed approach to be a valuable tool for the curation of genome-scale metabolic models. the outcome of its application to the genome-scale model b. cuenoti icg is a more accurate model version named as b. cuenoti imp 