BACKGROUND
in their natural environment, bacteria are confronted with a wide range of environmental conditions that change over time.
it is generally observed that, when facing a nutritionally rich environment, bacteria first use the “preferred” compound  and only later start metabolizing the other one.
the decrease in concentration of these compounds corresponds with the start of usage of the others.
in this respect, monod discovered the phenomenon of “diauxie”, i.e.
the microbial capability of utilizing the various nutrients regardless of their abundance but regulating their uptake through the modulation of the overall enzymatic state  <cit> .
this kind of switching of metabolic preference is characteristic of systems that optimize fitness  <cit> .
it is usually observed that, as the bacterium changes from one carbon source to another, growth is temporary halted, while a new of enzymes needed to metabolize alternative nutrients are synthesized.
this cellular regulation is likely to play a key role, i.e.
adjust cellular nutrient fluxes across the entire bacterial metabolic network as to produce the optimal growth rate  <cit> .
in other words, regardless of the c-source that is used in a specific growth phase, the metabolic network of an organism has to continuously and dynamically adjust to optimally sustain cellular growth.
accordingly, the optimal flux distribution achieved by growing cells sometimes changes discontinuously as the composition of the growth medium is varied.
therefore, a small change in nutrient concentration due, for example, to compounds exhaustion, may sometimes induce a large change in the enzymatic composition of the bacterium  <cit> .
examples of system-level cellular adjustments following changes in nutrients availability have been recently described exploiting –omics technologies  <cit> .
these works revealed presence of a deep and systemic re-organization of the overall microbial metabolic network in response to a variation in the composition/concentration of the surrounding nutrients.
the range and the entity of such modifications in organisms other than a few model microbes  are still unknown, despite they might provide interesting hints in the understanding of non-trivial metabolic phenotypes .
in order to systemically interpret such complex behaviours, scientists have started exploiting computational models  <cit>  and, as in  <cit> , their integration with phenotypic and gene expression data.
in this context, constraints-based approaches, in particular flux balance analysis , have been shown to be predictive of growth phenotypes  <cit>  and can be used to construct large scale metabolic models based on genome sequences and to infer the metabolic impact of perturbations in the external conditions   <cit> .
although the presence of metabolic switches in other organisms is less documented, a valuable exception is represented by the antarctic bacterium pseudoalteromonas haloplanktis tac <dig> .
phtac <dig> has been isolated from sea water sampled along the antarctic ice-shell and is considered one of the model organisms of cold-adapted bacteria.
it is capable of thriving in a wide temperature range and it has been suggested as an alternative host for the soluble overproduction of heterologous proteins, given its ability to rapidly multiply at low temperatures .
a genome-scale metabolic reconstruction of phtac <dig> has been recently published by our group  <cit> , allowing a system-level representation of its phenotypic landscape and boosting a holistic comprehension of its metabolic features.
interestingly, the growth curve of phtac <dig> in a complex medium  is characterised by a number of metabolic switches among the amino acids available in the medium  <cit> .
the progressive exhaustion of the different sets of amino acids  defines a hierarchical usage of carbon sources and, consequently, multiple phases in its growth curve, characterized by a different set of metabolized substrates.
as mentioned previously, this behaviour points towards a  reprogramming of cellular metabolism, the effects of which cannot be appreciated just looking at the resulting growth curve.
the aim of this study is to systemically investigate the metabolic adjustments that are  predicted to occur over time in a microbial cell grown in a complex medium, with particular attention to the global effects produced by multiple substrates switching .
constraint-based metabolic modelling was used to simulate microbial growth in a complex medium over several time steps, each characterized by a specific set of amino acids to be used as carbon and energy source.
as said, we examined these features in phtac <dig> taking advantage of its recently published genome-scale metabolic reconstruction and the availability of accurate physiological data in complex media.
results obtained from our simulations allow us to hypothesize the occurrence of a deep metabolic reprogramming following each transition in the availability of nutrients, with more than one half of its metabolic reactions involved.
interestingly, catabolic pathways are predicted not to be the only processes affected by such changes, with many central pathways that seem to be affected by apparent minor changes in the metabolized substrates.
furthermore, we show that our modelling framework is able to capture possible functional patterns and/or common regulatory features of the genes embedded in our reconstruction .
finally, we used our modelling framework to explore the effect of simulating a sub-optimal microbial growth, i.e.
accounting for the re-direction of part of the cellular resources in alternative metabolic processes at the expense of optimal growth rates.
the moma-based  <cit>  approach introduced here  was compared to fba-derived predictions, revealing interesting insights into the resulting predictions of fluxes distribution.
methods
model parameterization
the metabolic model used in this manuscript has been described in fondi et al.
<cit>  and can be found in additional file  <dig>  together with the codes used to produce the results reported herein.
the physiological data of phtac <dig> growth in soy-peptone complex medium  <cit>  was used to identify time steps corresponding to amino acid depletions.
overall, twelve different  time steps were identified along the complete growth curve and, for each of these time steps, the uptake rates of each amino acid were estimated based on their relative decrease in concentration and on the growth rate.
more in detail, to fit our model to wilmes et al.
data, first we computed the phtac <dig> yield.
to do this, we had to derive the difference in biomass for each time step.
first, we retrieved the optical density  and growth rate  at each time , deriving the ratio between biomass  and od, by dividing the reported weight at the 10 h mark  for the corresponding optical density  times the initial operation volume of the growth experiment .
then, we computed the difference  in biomass between each hour  and the following one .
we also computed, for each amino acid, the mass concentration  at each time step from the corresponding molar concentration values.
finally, we computed for each amino acid the difference in concentration  for each time step, considering this value as  <dig> when a concentration was increasing from an hour to the following one .
thus, for each time step, the biomass delta value represents the relative increase in biomass, while the sum of the amino acid delta values  represents the total concentration of carbon consumed in that time step.
the ratio of these values allowed us to compute the yield value corresponding to each time step.
to allow for growth simulations consistent with the experiment by wilmes et al.
<cit> , we also had to compute the uptake fluxes  for each amino acid at each time step.
to do this, we used the following formula: ufaai=concaai∑∀a∈aaconcai*griyieldi*1mwaa*10− <dig>
where: ufaai is the uf  of the amino acid aa  in the time step i; gri and yieldi are the phtac <dig> growth rate and yield, respectively, at the time step i; concaai and mwaa are the concentration  and molar weight values, respectively, of the amino acid aa in the time step i; ∑∀ a ∈ aaconcai is the summation of each amino acid concentration in the time step i. finally, we multiplied each uptake flux for a constant specific to the corresponding time step to allow for comparable flux values for each time step.
in silico growth was simulated for each phase using fba  with salts uptake as defined as in fondi et al.
<cit>  and adding amino acid uptake reactions with lower bound equal to the uptake rate estimated from growth data as described above.
modelling procedures
flux balance analysis  was employed to simulate fluxes distribution at each different time step  <cit> .
as during the last two time steps there was no growth and no amino acids uptake, we limited our analysis to the first ten time steps identified as previously described.
the reconstructed model was analysed using cobratoolbox- <dig> .
<dig>  <cit>  in matlab® r2012b .
gurobi  <dig>   solver was used for computational simulations presented herein.
the cobrapy package was also used during model expansion and preliminary growth simulations  <cit> .
statistical analyses on the predicted fluxes were performed using the r package  <cit> .
flux variability analysis  was performed using the cobra toolbox.
the ratio between maximum and minimum admissible flux of each reaction  in each of the growth phases was computed as: vi=fmax,i/fmin,i
with fmax,i and fmin.i representing the maximum and minimum admissible fluxes of the i
th reaction according to fva analysis, respectively.
on the basis of these spans, we determine the fixed  and flexible  parts of the metabolic network while achieving a particular metabolic objective .
identification of co-varying reactions and functional association among their genes
to identify co-varying reactions, we first removed reactions showing a constant trend throughout the simulation.
then, for each of the remaining reactions, we computed the difference  between the absolute values of the fluxes following each transition.
formally: dz=|fj,i+1|−|fj,i|
with fj,i representing the flux of the j
th reaction of the model in the i
th time step and dz the difference of such fluxes in the z
th transition.
as we divided the growth of phtac <dig> on complex medium in  <dig> different phases, for each reaction a vector embedding nine values of d was obtained.
pearson correlation was then computed for each pair of reactions.
finally, we extracted groups co-varying reactions by selecting those reactions sharing a pearson correlation value greater than  <dig> .
the sequences of these gene groups were queried to the string database  <cit>  to retrieve their possible functional association using its advanced programming inter-face .
string combined score includes in a single value  different hints of functional association  among those genes or among their orthologs in organisms other than the selected one .
nutritional-moma
with nutritional-moma we refer to the use of moma  <cit>  optimization to minimize the metabolic adjustments required at each  transition of the entire growth period analysed.
in its canonical formulation, the moma algorithm requires two models, the so called “wild-type” model and the “mutant” model.
moma can then be used to determine the flux distribution for the “mutant” model that minimizes the difference between the “mutant” model itself and the wild-type solution.
we have modified this formulation using, for each transition t among a given growth phase p and the previous one , the model of phase p as the “mutant” model and the model of phase p- <dig> as the wild type model.
the two models will different in the set of nutrients available during the computation of the flux distribution.
accordingly, this approach should provide a solution in which fluxes are computed in such a way that a minimal metabolic adjustment is computed when the model is presented to a different set of nutrients.
for sake of clarity, we named this alternative use of moma as nutritional-moma.
it is to be noticed that, as data from a t <dig> point  were not available, our nutritional-moma computations started from phase p <dig> .
as a result, plots for this set of analyses embed nine time points rather than ten as in fba simulations.
regulons and conserved motifs identification
the regprecise database  <cit>  was used to inspect the presence and the structure of putative regulons in phtac <dig>  specific upstream motifs finding was performed using the meme suite  <cit> .
genome-wide regulatory motifs searches were performed using the dminda server  <cit>  on the set of phtac <dig> upstream regions present in the door database  <cit> .
RESULTS
modelling procedures
we have recently reconstructed a genome-scale metabolic model of phtac <dig>  using it for inferring the metabolic adjustments of this bacterium induced by changes in gene expression following a temperature downgrade in this bacterium  <cit> .
here we used this model to investigate the metabolic rearrangements occurring during growth in a nutritionally rich medium.
a detailed analysis of the data reported by wilmes et al.
<cit>  allowed the identification of twelve distinct phases in the growth of phtac <dig> on peptone medium, each of them corresponding to a time step of one hour.
as during the last couple of phases almost no nutrients uptake was recorded, we limited our analyses to the first ten phases .
for each of these phases we identified the specific uptake rates of the different compounds  present in the growth medium  and/or  switches in the use of the available c-sources.
these values were used as input for ten different fba simulations  to derive the most likely fluxes distribution in the phtac <dig> metabolic model in each time step.
this allowed a system-level characterization of the metabolic changes occurring after variations in the usage  of the different carbon sources.
a schematic representation of the nutrients provided to the model in each of the time steps is reported in fig.
1a.fig.
<dig> summary of phtac <dig> genome-scale reprogramming following nutrients switching.
a. the nutrients provided to the model in each different growth phase according to  <cit>  b. heat map with log values of fluxes across all the phases.
c. number of flux carrying reactions in each growth phase.
d. number of flux-changing reactions in each growth phase.
the dashed line represents the average number of reactions carrying flux over all time points.
e. number of reactions whose flux is predicted to increase  and decrease  following each shift in the nutrients provided; the black line accounts for those reactions whose flux is predicted to decrease not as an effect of an imposed reduced growth rate during simulations.
f. functional annotation of reactions varying their flux across all the phases
prediction of core and switching reactions
first, we evaluated the predicted flux for each reaction in the model, across all the time points.
overall, we found that  <dig> reactions did not carry flux in our model in any of the growth phases; these reactions are probably not essential to sustain phtac125’s growth in the simulated conditions and were discarded from the following analyses.
conversely,  <dig> reactions were predicted to carry flux in at least one time point  and were considered for further analyses.
on average, excluding the last time point in which all the reactions are predicted to be turned off ,  <dig> reactions are predicted to carry flux for each of the analysed phases .
interestingly, this value is quite constant throughout the time points , revealing that, according to our simulations, the number of reactions necessary to support phtac <dig> growth is predicted to be somehow independent from the carbon sources  used as c and energy sources.
despite this apparently constant trend in the number of metabolic reactions employed by phtac <dig> in a complex medium, the shift among the different growth phases is presumably characterized by a relatively high number of reactions that are predicted  to change their flux .
on average,  <dig> reactions display at least a single variation in their predicted flux across all the time points .
however, in this case, a greater variability across all the time points is observed in our simulations .
gene-protein-reaction  rules of our model indicate that this latter set of reactions is encoded, on average, by  <dig> genes.
interestingly this number resembles that found when modelling the metabolic switch of s. coelicolor   <cit>  and represents around 12% of the phtac <dig> coding sequences .
according to our simulations, switching among the available carbon sources in a complex medium may have an impact on the overall metabolic network of phtac <dig> with more than half of flux-carrying reactions influenced by a change in the utilized c-source or in its uptake rate.
more in detail, we found that, on average,  <dig> reactions are predicted to decrease their flux throughout the growth period , whereas  <dig>  are predicted to increase it.
however, the number of reactions whose flux decreases might be biased due to a general decrease, across the time steps, of the growth rate.
more specifically, this systemic bias is related to the amino acid uptake rates, which have been derived from the physiological data by wilmes et al.
<cit> , whose constant decrease lead, for some reactions, to a flux reduction in each time step.
thus, we adjusted the set of “decreasing reactions”, by removing those for which we observed a consistent decrease for all the growth phases.
although the normalization did not affect the general trend, in that the normalized set  and the not normalized one  have similar trends, the number of decreasing reactions is  comparable to that of the increasing reactions.
furthermore, we also computed a normalized flux distribution for each of the modelled growth phases, expressing them as a fraction of the predicted growth rate and evaluated whether this procedure led to different results .
results  revealed no major differences in the overall trends compared to the original calculation of fluxes distribution, suggesting that our results hold true regardless of the normalization procedure adopted to account for the different  growth rates predicted by the model across the growth period.
not all phtac <dig> metabolic pathways are impacted by this switching of nutrients in the same manner, according to our simulations.
figure 2a shows the hypothetical number of flux-carrying reactions for five pathways, i.e.
tca cycle, lys biosynthesis, glu metabolism, val, leu and ile biosynthesis and degradation .
tca cycle, for example, displays and increasing trend in the number of flux-carrying reactions according to our modelling framework; this is consistent with the exhaustion of amino acids  whose degradation provides important tca cycle intermediates, i.e.
oxaloacetate, fumarate and 2-oxo-glutarate and, consequently, with the necessity to activate those reactions leading to the biosynthesis of such compounds.
conversely, lys and glu metabolic routes display an overall constant trend , with a similar number of active reactions across the different simulated growth phases.
this is in line with i) the necessity to use  of the lysine biosynthetic route to synthesize diaminopimelic acid, an essential component of bacterial cell wall  and ii) with the importance of glu metabolism for phtac <dig> .
finally, val, leu and ile biosynthesis and degradation pathways display an opposite trend one another .
intuitively, this reflects the necessity to synthesize these molecules in the first part of the growth phase  and the necessity to catabolise them once phtac <dig> is using those amino acids as carbon sources, respectively.fig.
<dig> changes in the central metabolism of phtac <dig>  a. the number of active  reactions for five major pathways across all the time points is shown.
b. a simplified representation of the interconnections in the central metabolism of phtac <dig>  dashed lines indicate the presence of more than one reaction between the connected compounds.
modified from  <cit>
we next evaluated whether changes in the kind and number of utilized substrates across all the time points was also parallel to major structural changes at the whole metabolic network level.
we used flux variability analysis  to perform a comprehensive exploration of alternate optimal routes in phtac <dig> metabolic network, in each of the predicted growth phases .
this analysis revealed that the number of “flexible reactions” is predicted to remain constant throughout all the growth phases , suggesting that the plasticity of phtac <dig> metabolic network is not particularly dependent on the set of metabolized substrates.
this suggests the presence of an efficient adaptation of phtac125’s metabolic network to face a nutritionally rich but probably variable environment.
overall, according to our model, the central phases of the growth curve are characterized by the largest set of reactions changing their flux , consistently with the major changes in the set of utilized nutrients observed experimentally by wilmes et al.
<cit>  .
the first important switch in the panel of utilized compounds occurs between phase  <dig> and  <dig> .
according to experimental evidence from wilmes et al.
<cit> , in this time point phtac <dig> stops using ser, asn and asp as c and n sources, relying only on glutamate for sustaining growth.
the importance of glutamate was previously observed from fed-batch cultivation experiments as it was the most strongly metabolised amino acid in all growth experiments  <cit> .
also, glu  is the amino acid allowing the fastest growth rate  among all the  <dig> amino acids when growth in twenty minimal growth media  * h) was simulated .
during p <dig>  phtac <dig> relies only on glutamate for sustaining growth and needs to utilize this compound to derive its main building blocks.
l-glutamate:nadp+ oxidoreductas  is predicted to be a key reaction in this phase as it allows the conversion of glutamate to 2-oxoglutarate to be used as energy source feeding the tca cycle .
more in general, our modelling framework predicts a deep re-programming at the whole metabolic network level following this shift, with  <dig> reactions changing their predicted flux; more in detail,  <dig> reactions are predicted to decrease their flux , whereas  <dig> are predicted to carry an increased flux .
the remaining  reactions are predicted to operate in the reverse direction compared to the previous growth phase following the switch.
among the first set of reactions,  <dig> are predicted to be turned on in respect to the previous phase, whereas  <dig> are predicted to be turned completely off.
this latter set of reactions includes, as it might be expected, the transporters of serine, asparagine and aspartate, i.e.
those compounds that are no longer available in the medium and that cannot be internalized anymore.
the cell faces the absence of some of these compounds by redirecting  of the glutamate available in the medium to their synthesis.
indeed, for example, the reaction encoded by l-aspartate:l-glutamine amido-ligase  leading to the formation of asparagine from gln  is predicted to double its flux after this metabolic transition .
the exhaustion of asp is probably compensated by reversing reaction r <dig> , generating asp from glu .
furthermore, our simulation indicated a drop to no flux in the reaction encoded by threonine dehydratase  and catalysing the conversion of thr to production of 2-oxobutanoate , a precursor of val, leu and ile .
a decrease in the availability of asp  might impair the biosynthesis of 2-oxobutanoate following this pathway.
our model predicts that phtac <dig> faces this perturbation by reversing the flux of reaction r <dig> , leading to the production of 2-oxobutanoate  from o-succinyl-l-homoserine.
after the following switch , wilmes et al.
observed that phtac <dig> starts using ala, leu, gly and tyr  .
according to our modelling framework, this involves a predicted change in the flux of  <dig> reactions , with  <dig>   <dig> and  <dig> of them decreasing, increasing or changing the direction of their flux, respectively.
consistently with this new set of amino acids, the model predicts that most of the biosynthetic routes leading to the production of such compounds result turned “off”.
this is the case of the path leading to the synthesis of leu in the branched-chain amino acid biosynthetic route .
the same occurs for ala biosynthesis  and l-aspartate 4-carboxy-lyase  are predicted to be turned off) and lys biosyntheses.
interestingly, in the latter case, only the final step of the whole pathway  is predicted to carry no flux, consistently with the observation that the rest of the pathway is required for the synthesis of meso- <dig> -diaminopimelate, an essential precursor for peptidoglycan assembly .
amino acid biosynthetic routes are not the only pathways affected by this metabolic switch, according to our simulated growth model.
two tca cycle reactions are predicted to be turned on following the utilization of ala, leu, gly and tyr, i.e.
those leading to the production of 2-oxoglutarate from isocitrate .
this raises the intriguing question on the source of 2-oxoglutarate in the previous growth phases given that those tca reactions were predicted to carry no flux.
a possible explanation is provided by the observation that 2-oxoglutarate can be obtained from the carbon skeletons of several five-carbon amino acids through a first conversion into glu, which is then oxidatively deaminated by glutamate dehydrogenase to yield α-ketoglutarate.
however, the presence of glu in the medium would not favour this solution as the major fraction of glu can be obtained without additional energy expenses.
alternatively, 2-oxoglutarate can be obtained from the conversion of glu  to tyr .
in our model this reaction  is predicted to be “on” during p <dig> thus allowing the synthesis of 2-oxoglutarate .
however, starting to use tyr present in the medium during p <dig> might cause reaction r <dig> to be turned off and, consequently, the necessity to synthesize 2-oxoglutarate from succinate to maintain the functioning of tca cycle.
similarly, the production of fumarate from malate  is predicted to carry no flux in the shift to p <dig> .
again, this is consistent with the degradation of amino acids as a major source of important metabolic intermediates, as fumarate can be obtained from the degradation of aromatic amino acids as tyr, available in the medium during p <dig> and then tunnelled into the tca cycle.
the time-resolved growth data from wilmes et al.
<cit>  show that, after other minor transition in which few additional amino acids are degraded , the final switch involved the utilization of histidine as the only c source and causes a major drop in the growth rate of phtac <dig>  <cit> .
in our model, this corresponds to the highest number of flux-changing reactions  and thus, to the deepest reprogramming encountered by phtac <dig> during this growth curve.
as it might be expected, histidine biosynthetic reactions  are predicted to carry zero flux following this transition.
despite the entire pathway is predicted to carry no flux, aicar -5-amino-4-imidazolecarboxamide), one of the intermediates of histidine biosynthesis and a crucial precursor in purine metabolism, might still be synthesized through reaction r <dig> .
this compound, thanks to the flux predicted to be carried in this phase by reaction r <dig>  leads to the synthesis of 5-phospho-alpha-d-ribose 1-diphosphate from d-ribose 5-phosphate.
furthermore, the switch to the utilization of his as the sole c-source is predicted to be responsible for the stop of the catabolic routes of the previously degraded amino acids  and, consequently, in the stop of the production of important cellular intermediates .
according to our simulation, this causes the re-activation of key metabolic reactions that cannot rely on many degradation intermediates as in the previous phase.
these include those belonging to the tca cycle  and purine metabolism .
finally, the reactions involved in the conversion  of his into glu  are predicted to be turned on following this final transition, allowing the production of glutamate and, from these, all the major pathways found to carry flux also in phase p <dig> .
a sub-optimal objective function predicts alternative fluxes distribution
it is worth noticing that all the simulations described up to now were conducted under the fba canonical assumption of biomass optimality, i.e.
assuming that all metabolic fluxes in the cell are geared towards the production of biomass in each moment of the growth curve.
however, situations in which cells invest substantial resources in alternative metabolic processes at the expense of optimal growth rates have been analysed and described quite extensively.
according to this scenario, cells may invest substantial  resources in a specific metabolic process at the expense of optimal growth, this being reflected by sub-optimal flux distributions  <cit> .
this situation might be observed, for example, when microbes allocate energetic resources for anticipating changing environmental conditions at the expense of optimal growth  <cit> .
also, when cells are exposed to nutrients fluctuations they might respond with a minimal metabolic adjustment, to avoid the waste of protein synthesis and degradation necessary to reprogram the entire metabolic network and to simultaneously achieve two objectives, i.e.
rapid and minimal adjustments.
this latter scenario may indeed resemble the actual competition for nutrients that emerges in natural environments and that has been proposed to be reflected by the sequential uptake of nutrients  <cit> .
to explore this alternative scenario, we have accounted for the possibility that the phtac <dig> biomass objective function could not be fully optimized but, instead, in a near-optimal or sub-optimal state.
since fba classically assumes the optimization of the biomass production flux to compute the most likely fluxes distribution inside the cell, it cannot per se provide hints concerning alternative  fluxes distribution.
for this reason, we have used moma  <cit>  to study the hypothetical fluxes distribution when minimizing the metabolic adjustments required at each  transition of the entire growth period analysed .
this differs from the canonical formulation of moma in which the effects of a gene knock-out are evaluated by providing an approximate solution for a sub-optimal growth flux state , nearest in flux distribution to the unperturbed state .
applying this modelling strategy to our study case revealed that the choice of the optimization criterion  has a great influence on the predicted fluxes distribution, both in terms of the fraction of the entire network required to sustain phtac <dig> growth and on the set of active pathways.
first, a higher number of active reactions compared to the original fba predictions were predicted for each of the growth phases, when the implemented nutritional-moma approach was used .
this gap is even more evident in the first seven growth phases, whereas the difference between the two approaches becomes negligible in the last three phases.
however, despite the two approaches predict a similar number of active reactions in these two final time points, the two sets of reactions appear to be quite different .fig.
<dig> comparison between nutritional-moma and fba.
here we show the number of predicted flux carrying reactions in each growth phase for fba  and nutritional-moma  optimization on the phtac <dig> model.
also, the number of shared reactions identified by the two approaches is shown
the fact that nutritional-moma always predicts a higher number of active reactions might be explained by considering that to optimize the model using fba identifies the flux distribution maximizing the objective function , regardless of the metabolic states in the previous time steps.
conversely, nutritional-moma will seek for the solution  that is closer to the one of each previous time step .
therefore, the results of moma will likely include a higher number of reactions in the model to produce biomass, since it will minimize the changes  with respect to the previous time steps, while activating novel reactions to cope with the changes in the nutrient composition.
we then investigated whether these differences in the sets of predicted flux-carrying  reactions between the two approaches were involved in specific pathways or, rather, spanned over a larger part of the entire phtac <dig> metabolic network.
thus, for each metabolic pathway in the model, we computed the fraction of reactions predicted to be active by nutritional-moma, fba and both methods .
results of this analysis revealed that the choice of the optimization method impacts the predicted fluxes distribution not only for what concerns the activity of the peripheral  pathways, i.e.
those pathways that start the degradation of amino acids and then tunnel them into the central metabolism.
indeed, we found many central processes  in which the proportion of active reactions is  specific for each of the two optimization strategy .
notably, for the first six growth phases, most of the reactions predicted to be active under fba optimizations  are carrying flux also adopting the nutritional-moma approach.
in some cases, however, entire pathways are differently predicted to be active/inactive by the two approaches .
the last two phases of the growth nutritional-moma and fba predict a similar number of active reactions but, from a functional viewpoint, deep differences exist in that entire pathways are predicted to be active only under a specific optimization method .fig.
<dig> functional differences between nutritional-moma and fba predictions.
here we show the proportion of reactions predicted to be active by nutritional-moma , fba  and both methods  for each main functional category represented in the phtac <dig> reconstruction
to summarize, using our moma-based approach we identified a higher number of reactions predicted to be carrying flux in respect to the fba-based optimization.
the nutritional-moma approach predicted larger active metabolic networks throughout the analysed phases possessing, on average, about  <dig> flux-carrying reactions more than the one simulated with fba.
furthermore, despite most of the active reactions identified through fba were also identified by the nutritional-moma approach, in some cases these predictions differed significantly from a functional viewpoint .
at present, further experimental evidences are needed to shed light on the real number and function of reactions used by phtac <dig> in each time point  and, in other words, to infer how far from the actual fluxes distribution our in silico predictions are.fig.
<dig> flux correlation analysis.
heatmap accounting for the pearson correlation of all the flux difference vectors across all the time points.
the metabolic process of each reaction is also reported
flux correlation analysis identifies functionally associated genes
in this section we analyse more in depth the flux matrix .
this matrix contains the fba-predicted flux of the reactions  for each of the time steps .
as such, it allows capturing the co-variation of reactions in the model; in other words, we identified as paired reactions those whose predicted fluxes displayed a similar trend throughout the time points.
these, in turn, may represent functional partners in the cell and/or may belong to the same metabolic pathway/module.
we represented co-varying reactions  in the form of a heatmap as shown in fig.
<dig>  in most cases, the small clusters of reactions showing high pearson correlation values are involved in the same metabolic process.
this is compatible with a scenario in which, following a nutrient switch, entire metabolic pathways are activated  to face the novel environmental condition .
nevertheless, grouping reactions with a lower but still significant threshold , clusters start to embed reactions belonging to other metabolic processes.
these, in turn, may represent previously undetected functional associations between genes and/or entire pathways.
overall, our method led to the identified  <dig> different clusters fig.
6a, comprising  <dig> reactions.
according to the gpr of our model, these reactions were encoded by  <dig> genes.fig.
<dig>
a. co-varying reactions clusters identification.
common flux trends  for the reactions embedded in each of the  <dig> clusters.
b. the distribution of string combined scores among all the genes embedded in each cluster of genes  and the number of genes embedded by each cluster .
the grey line represents the median of the combined score computed for each possible pair of genes in the model
we analysed these clusters to answer two biological questions, namely whether the embedded genes i) are de facto functionally related and ii) share conserved upstream motifs .
groups of co-varying genes display higher combined evidence scores
to address the first point, we exploited the combined evidence score provided by the string database.
as shown in fig.
6b, the median of the combined score for the genes embedded in the same cluster is, in most cases, well above the median of all the possible combinations of the genes embedded in the model, suggesting that flux coupling in the model occurs among bona fide functional patterns.
accordingly, this result corroborates the capability of our modelling framework to identify functional metabolic modules and the functional associations of their encoding genes in phtac <dig>  interestingly, in some cases apparently un-related genes  seem to be embedded in the clusters.
despite these instances may represent erroneous predictions of our model, they might also suggest the presence of still undetected functional associations and/or the use of common pools of chemical intermediates by the corresponding reactions .
we here present a description of the functional relationships retrieved for genes belonging embedded in three of these clusters .
a detailed list of the reactions and genes embedded in all the clusters is provided as additional file  <dig> to allow further experimentation and/or in silico analyses.fig.
<dig> sample string clusters.
evidence network, co-expression instances and co-occurrence patterns for clusters  <dig> ,  <dig>  and  <dig> .
red asterisks in g indicate those genes known to belong to the same regulon according to regprecise database
cluster  <dig>  this cluster is composed of seven reactions encoded by  <dig> genes.
three of these reactions are predicted to be catalysed by enzymes encoded by multiple genes, i.e.
cytochrome-c oxidase encoded by pshaa2869- <dig>  ferrocytochrome-coxygen oxidoreductase encoded by pshaa1842- <dig> and atp synthase encoded by pshaa3007- <dig>  according to the gpr rules of our model, the simultaneous presence of all the corresponding genes is required to carry flux of each of these reactions.
the functional association network obtained probing the string database with all the genes of cluster  <dig> is reported in fig.
7a.
genes encoding proteins involved in the same enzymatic reaction appear to be highly interconnected among themselves with evidences that included co-expression, gene fusions, and experimental evidences .
this is not surprising, as it might be expected that genes encoding proteins that participate in the same enzymatic complex should display coupled regulation and co-occurrence patterns in other organisms.
links, however, are also found among members of different enzymatic complexes, consistent with the fact that these reactions belong to the overall process of oxidative phosphorylation.
the ndk gene  is also found connected to one of the sub-cluster fig.
7a).
this protein plays a major role in the synthesis of nucleoside triphosphates and is involved in purine and pyrimidine metabolism.
the coupling of this reaction’s flux with those belonging to oxidative phosphorylation can thus be explained by their common final products , as optimizing the model for biomass production would probably force the flux through both these pathways.
finally, the fact that orthologs of these genes are found to be co-expressed and located in physical proximity in other  genomes  support the functional association predicted by our model.
cluster  <dig>  we found eight reactions  embedded in this cluster , involved in the following metabolic pathways: glycerolipid metabolism, pentose phosphate pathway and glycolysis.
half of the reactions of this cluster are involved in the conversion of glycerol into glycerate 3-pshosphate that, in turn, is connected to glycolysis/gluconeogenesis.
as shown by the string analysis output, many evidences exist of a functional association among the genes of this cluster, in most cases including instances of co-expression and physical proximity fig.
7e and f, respectively).
two genes of the cluster remained disconnected in the string network , both encoding a glycerol kinase and assigned to the same reaction in our metabolic reconstruction .
accordingly, despite no previously detected functional associations were retrieved between these two genes and the other embedded in the same cluster, their presence could be accounted for by the use of common precursors  with the reactions encoded by cluster  <dig> genes.
indeed, fluctuations in the availability of such compound  during the simulated growth transitions explains the coupling of such reactions.
cluster  <dig>  this cluster embeds nine reactions encoded by  <dig> genes involved in arg biosynthetic pathway argah, b, c, d, e and in the metabolism of pro ).
the string-based evidence network in fig.
7g revealed hints on the functional associations existing among the genes of the cluster including co-occurrence of orthologs of these genes in other  genomes, as well as co-expression and conserved genomic neighbourhoods fig.
7h and i, respectively).
the genes embedded in this cluster encode enzymes involved in the the formation of arg and pro from glu, using asp and then orn as intermediates.
these reactions appear to carry flux especially in the first three phases of the growth , that is when aspartate  is present in the medium.
afterwards, a decrease in the flux carried by these reactions is observed , followed by an overall constant flux trend in the next transitions.
accordingly, the correlation among the fluxes of these reactions could be explained by the need to turn these pathways “on”  when glutamate and required intermediates are present  in the growth medium, respectively.
a gene encoding an acetyl-coa synthetase  is also embedded in this cluster, despite it has not connections with arginine and proline metabolism-related genes.
however, its presence might reflect the necessity to synthesize acetyl-coa from acetate during the first four growth phases identified, i.e.
in the absence of ketogenic amino acids.
once phtac <dig> starts metabolising the second set of amino acids, acetyl-coa can be obtained from the degradation of specific amino acids  and thus the reaction encoded by acsa  is predicted to stop carrying flux.
analysis of upstream sequences identifies putative regulation mechanisms
in the previous section we reported a list of gene sets which are, presumably, functional partners and whose activity is concerted during the different metabolic switches.
if this holds, it is reasonable to hypothesize that some of these clusters might embed co-regulated genes that, in turn, may share conserved upstream motifs .
thus, to address this point we applied an ad hoc computational pipeline  to analyse the upstream regions of the genes  embedded in each cluster.
this resulted in the identification of motifs that are shared in the upstream regions of all the genes for each cluster.
upstream regions were identified using door   <cit> .
these have been annotated using three different motif databases , to identify motifs with putatively related functions and that might be recognized by putatively related transcription factors .
the whole set of annotations has been manually curated in a conservative way to obtain high-confidence motif annotations.
finally, the obtained putative tfs have been cross-validated by comparing their associated functions with those related to the genes embedded in the clusters.
the whole set of data related to the identified motifs and their annotations , for each database, is provided as additional file  <dig>  on average, we identified for each cluster a large number of annotations according to the different databases .
however, after manual curation, the number of annotations was reduced to  <dig>  divided in  <dig> clusters .
most of the remaining annotations  come from prodoric, while the other databases have a minor contribution of high-quality annotations .
for each cluster, we eventually compared the putative tfs with the functions of the genes embedded in the cluster, using both literature information and the regprecise database.
with few exceptions that will be described in details, we found no correlation between the putative tfs and the regulation/biological function associated to the genes embedded in the clusters, suggesting that the majority of the co-varying groups share a common regulation but are regulated by a number of different tfs.
we are describing below the clusters for which we found agreement between the predicted tfs and the function of the embedded genes .table  <dig> main features of the putatively co-regulated clusters found during flux-correlation analysis.
in this table, for each cluster, we report its number, the name of the regulator identified, the genes embedded in it and the conserved motif found upstream of its genes
cluster  <dig>  this cluster, as described previously, is involved in the metabolism of arg and pro.
a search on the regprecise regulon database  <cit>  revealed that, at least in the case of arg, the genes embedded in the cluster of fig.
7g are part of a common regulon .
the only exception is represented by argd which, according to the regprecise database is not part of the arginine regulon in phtac <dig>  and the other genes involved in the biosynthesis of pro.
however, both the computational approach described above and a manual search for putative transcription factor binding sites upstream of these genes detected the already known conserved motif shared by the other arg genes of the regulon .
this suggests that the all the genes are part of the arginine regulon in phtac <dig>  more in general, this finding confirms that, starting from our modelling outcomes it is possible to infer biologically consistent patterns of gene co-expression.
cluster  <dig>  this cluster includes genes from the sdh operon , encoding the succinate dehydrogenases system, as well as other enzymes involved in energy conversion, such as other dehydrogenases  and other enzymes with a role in nucleotide modifications .
our computational pipeline predicted two different regulator genes associated with this cluster, i.e.
ccpa and psra .
these regulators encode proteins are involved in the carbon and fatty acid metabolism, respectively.
comparing the functional role of the tfs with the cluster annotation, we can easily relate the set of genes encoding dehydrogenases with ccpa and psra, which is also confirmed by regprecise database.
the other three genes are indirectly involved in the oxidative metabolism, in that they are involved in the metabolism of fad.
cluster  <dig>  this cluster includes five genes, with three of them involved in the synthesis of folate , while the remaining two encode a galactokinase  and a galactose-1-epimerase .
all these genes are involved in the metabolism of nucleotides, since folate is necessary to synthesize tdp sugars from udp, while the other genes encode enzymes that are involved in the synthesis of udp.
the analysis of the upstream sequence of these genes revealed a conserved motif similar to the galr binding site, thus suggesting that these genes are being co-regulated .
CONCLUSIONS
here we have used constraint-based metabolic modelling to provide a time-resolved, systems-level scheme of phtac <dig> metabolic re-programming following nutrients switching in a nutritionally complex medium.
such features have been analysed using the metabolic model and growth data of the antarctic bacteria p. haloplanktis tac <dig>  previous experimental tests revealed a number of nutrients switches in this microorganism when grown in a complex medium  <cit> , consistent with the fact that this antarctic organism is adapted to fast growth in a fairly rich  environment   <cit> .
indeed, sequential uptake of nutrients is thought to emerge when competition for nutrients is present  <cit> .
our modelling framework identified the central phase of the growth curve as a probable key reprogramming point, with more than  <dig> reactions predicted to adjust their flux in these time points.
this corresponds  to the exhaustion of most of the first metabolized c sources , a time step in which only glu is used as the sole c source, and the final part of the growth curve in which phtac <dig> stops using glu and relies on a completely different set of nutrients.
our model highlights the occurrence of such an adaption and the need for reprogramming a large set of reactions to maintain an efficient metabolic network.
a similar scenario is observed at the end of the growth, when his is used as the sole c source.
this transition is the most demanding, requiring a change in the predicted flux of almost  <dig> reactions.
taken together, our simulation indicates the presence of an almost constant number of reactions  that are required to sustain life across all the time points examined.
interestingly, this set of reactions resembles, in percentage, the size of the minimal metabolic network predicted to be active in e. coli  <cit>  .
to maintain such a constant trend, however, a deep reprogramming of the whole cellular metabolism appears to be necessary during the entire growth period, as up to more than  <dig> biochemical reaction, display at least one change in the predicted carried flux.
according to our simulations, these changes do not only involve peripheral metabolic pathways  but also a number of central pathways, e.g.
tca cycle, glycolysis and pp pathway.
tca cycle, for example, displays an increasing trend in the number of flux-carrying reactions parallel to the exhaustion  of the amino acids whose degradation provide key tca cycle intermediates.
this is in line with the assumption that fluxes distribution within the cell are influenced by the entry-point of the a given c-source into the metabolic network  <cit> .
indeed, the c sources provided in our simulations have different distances from the phtac <dig> central metabolism  and thus are expected to cause a re-wiring of an important fraction of the network.
this, in turn, suggests the presence of an efficient metabolic reprogramming machinery  to continuously and promptly adapt to this nutritionally changing environment and/or to the exhaustion of the preferred carbon source.
modelling and dividing the growth curve in discrete time points has allowed us to infer common trends in the predicted flux patterns of the reactions of the model.
we have shown that such coupled patterns likely correspond to reactions and enzymes that work in a concerted fashion and that, in some cases, represent functional partners and/or be encoded by co-regulated genes.
the modelling framework we have set up here has allowed gaining insights on the regulation of phtac <dig> metabolic genes and, as in the case of the arg regulon, to expand the current knowledge on commonly regulated genes.
furthermore, we have assembled a dataset of putatively co-regulated genes that can be used for further manual curation and/or subject to experimental validation.
despite fba-based predictions led to the identification of biologically consistent trends, the nutritional-moma approach we have implemented here suggests that care needs to be taken when choosing the objective function for constraint-based metabolic modelling.
indeed, when we accounted for the possibility that cells could be in a near-optimal or sub-optimal state  the two modelling frameworks predicted quite a different topological and functional reprogramming of phtac <dig> cells.
as it is currently not possible to unambiguously discern among these two alternative solutions provided by fba and nutritional-moma, further experimental evidence  of phtac <dig> cells grown in complex medium will allow deriving a clearer picture of those pathways that are really active during the growth and, consequently, which of the two approaches outperforms the other.
additional files
additional file 1: the metabolic reconstruction and the matlab codes used in this work.
additional file 2: additional analyses on the phtac <dig> metabolic reconstruction showing i) the effect of expressing fluxes as a fraction of the biomass production flux and ii) number of active reactions in each phtac <dig> metabolic pathway following the nutrients switching.
additional file 3: the list of the reactions and genes embedded in all the clusters identified through flux correlation analysis.
additional file 4: the whole set of data related to the identified motifs and their annotations , for each probed database.
abbreviations
fbaflux balance analysis
fvaflux variability analysis
momaminimization of metabolic adjustment
phtac125
pseudoalteromonas haloplanktis tac125
tcatricarboxylic acid
tftranscription factor

 in their natural environment, bacteria are confronted with a wide range of environmental conditions that change over time. it is generally observed that, when facing a nutritionally rich environment, bacteria first use the “preferred” compound and only later start metabolizing the other one. the decrease in concentration of these compounds corresponds with the start of usage of the others. in this respect, monod discovered the phenomenon of “diauxie”, i.e. the microbial capability of utilizing the various nutrients regardless of their abundance but regulating their uptake through the modulation of the overall enzymatic state. this kind of switching of metabolic preference is characteristic of systems that optimize fitness. it is usually observed that, as the bacterium changes from one carbon source to another, growth is temporary halted, while a new of enzymes needed to metabolize alternative nutrients are synthesized. this cellular regulation is likely to play a key role, i.e. adjust cellular nutrient fluxes across the entire bacterial metabolic network as to produce the optimal growth rate. in other words, regardless of the c-source that is used in a specific growth phase, the metabolic network of an organism has to continuously and dynamically adjust to optimally sustain cellular growth. accordingly, the optimal flux distribution achieved by growing cells sometimes changes discontinuously as the composition of the growth medium is varied. therefore, a small change in nutrient concentration due, for example, to compounds exhaustion, may sometimes induce a large change in the enzymatic composition of the bacterium. examples of system-level cellular adjustments following changes in nutrients availability have been recently described exploiting –omics technologies. these works revealed presence of a deep and systemic re-organization of the overall microbial metabolic network in response to a variation in the composition/concentration of the surrounding nutrients. the range and the entity of such modifications in organisms other than a few model microbes are still unknown, despite they might provide interesting hints in the understanding of non-trivial metabolic phenotypes. in order to systemically interpret such complex behaviours, scientists have started exploiting computational models and, as in, their integration with phenotypic and gene expression data. in this context, constraints-based approaches, in particular flux balance analysis, have been shown to be predictive of growth phenotypes and can be used to construct large scale metabolic models based on genome sequences and to infer the metabolic impact of perturbations in the external conditions. although the presence of metabolic switches in other organisms is less documented, a valuable exception is represented by the antarctic bacterium pseudoalteromonas haloplanktis tac. phtac has been isolated from sea water sampled along the antarctic ice-shell and is considered one of the model organisms of cold-adapted bacteria. it is capable of thriving in a wide temperature range and it has been suggested as an alternative host for the soluble overproduction of heterologous proteins, given its ability to rapidly multiply at low temperatures. a genome-scale metabolic reconstruction of phtac has been recently published by our group, allowing a system-level representation of its phenotypic landscape and boosting a holistic comprehension of its metabolic features. interestingly, the growth curve of phtac in a complex medium is characterised by a number of metabolic switches among the amino acids available in the medium. the progressive exhaustion of the different sets of amino acids defines a hierarchical usage of carbon sources and, consequently, multiple phases in its growth curve, characterized by a different set of metabolized substrates. as mentioned previously, this behaviour points towards a reprogramming of cellular metabolism, the effects of which cannot be appreciated just looking at the resulting growth curve. the aim of this study is to systemically investigate the metabolic adjustments that are predicted to occur over time in a microbial cell grown in a complex medium, with particular attention to the global effects produced by multiple substrates switching. constraint-based metabolic modelling was used to simulate microbial growth in a complex medium over several time steps, each characterized by a specific set of amino acids to be used as carbon and energy source. as said, we examined these features in phtac taking advantage of its recently published genome-scale metabolic reconstruction and the availability of accurate physiological data in complex media. results obtained from our simulations allow us to hypothesize the occurrence of a deep metabolic reprogramming following each transition in the availability of nutrients, with more than one half of its metabolic reactions involved. interestingly, catabolic pathways are predicted not to be the only processes affected by such changes, with many central pathways that seem to be affected by apparent minor changes in the metabolized substrates. furthermore, we show that our modelling framework is able to capture possible functional patterns and/or common regulatory features of the genes embedded in our reconstruction. finally, we used our modelling framework to explore the effect of simulating a sub-optimal microbial growth, i.e. accounting for the re-direction of part of the cellular resources in alternative metabolic processes at the expense of optimal growth rates. the moma-based approach introduced here was compared to fba-derived predictions, revealing interesting insights into the resulting predictions of fluxes distribution. methods model parameterization the metabolic model used in this manuscript has been described in fondi et al. and can be found in additional file together with the codes used to produce the results reported herein. the physiological data of phtac growth in soy-peptone complex medium was used to identify time steps corresponding to amino acid depletions. overall, twelve different time steps were identified along the complete growth curve and, for each of these time steps, the uptake rates of each amino acid were estimated based on their relative decrease in concentration and on the growth rate. more in detail, to fit our model to wilmes et al. data, first we computed the phtac yield. to do this, we had to derive the difference in biomass for each time step. first, we retrieved the optical density and growth rate at each time, deriving the ratio between biomass and od, by dividing the reported weight at the 10 h mark for the corresponding optical density times the initial operation volume of the growth experiment. then, we computed the difference in biomass between each hour and the following one. we also computed, for each amino acid, the mass concentration at each time step from the corresponding molar concentration values. finally, we computed for each amino acid the difference in concentration for each time step, considering this value as when a concentration was increasing from an hour to the following one. thus, for each time step, the biomass delta value represents the relative increase in biomass, while the sum of the amino acid delta values represents the total concentration of carbon consumed in that time step. the ratio of these values allowed us to compute the yield value corresponding to each time step. to allow for growth simulations consistent with the experiment by wilmes et al., we also had to compute the uptake fluxes for each amino acid at each time step. to do this, we used the following formula: ufaai=concaai∑∀a∈aaconcai*griyieldi*1mwaa*10− where: ufaai is the uf of the amino acid aa in the time step i; gri and yieldi are the phtac growth rate and yield, respectively, at the time step i; concaai and mwaa are the concentration and molar weight values, respectively, of the amino acid aa in the time step i; ∑∀ a ∈ aaconcai is the summation of each amino acid concentration in the time step i. finally, we multiplied each uptake flux for a constant specific to the corresponding time step to allow for comparable flux values for each time step. in silico growth was simulated for each phase using fba with salts uptake as defined as in fondi et al. and adding amino acid uptake reactions with lower bound equal to the uptake rate estimated from growth data as described above. modelling procedures flux balance analysis was employed to simulate fluxes distribution at each different time step. as during the last two time steps there was no growth and no amino acids uptake, we limited our analysis to the first ten time steps identified as previously described. the reconstructed model was analysed using cobratoolbox-. in matlab® r2012b. gurobi solver was used for computational simulations presented herein. the cobrapy package was also used during model expansion and preliminary growth simulations. statistical analyses on the predicted fluxes were performed using the r package. flux variability analysis was performed using the cobra toolbox. the ratio between maximum and minimum admissible flux of each reaction in each of the growth phases was computed as: vi=fmax,i/fmin,i with fmax,i and fmin.i representing the maximum and minimum admissible fluxes of the i th reaction according to fva analysis, respectively. on the basis of these spans, we determine the fixed and flexible parts of the metabolic network while achieving a particular metabolic objective. identification of co-varying reactions and functional association among their genes to identify co-varying reactions, we first removed reactions showing a constant trend throughout the simulation. then, for each of the remaining reactions, we computed the difference between the absolute values of the fluxes following each transition. formally: dz=|fj,i+1|−|fj,i| with fj,i representing the flux of the j th reaction of the model in the i th time step and dz the difference of such fluxes in the z th transition. as we divided the growth of phtac on complex medium in different phases, for each reaction a vector embedding nine values of d was obtained. pearson correlation was then computed for each pair of reactions. finally, we extracted groups co-varying reactions by selecting those reactions sharing a pearson correlation value greater than. the sequences of these gene groups were queried to the string database to retrieve their possible functional association using its advanced programming inter-face. string combined score includes in a single value different hints of functional association among those genes or among their orthologs in organisms other than the selected one. nutritional-moma with nutritional-moma we refer to the use of moma optimization to minimize the metabolic adjustments required at each transition of the entire growth period analysed. in its canonical formulation, the moma algorithm requires two models, the so called “wild-type” model and the “mutant” model. moma can then be used to determine the flux distribution for the “mutant” model that minimizes the difference between the “mutant” model itself and the wild-type solution. we have modified this formulation using, for each transition t among a given growth phase p and the previous one, the model of phase p as the “mutant” model and the model of phase p- as the wild type model. the two models will different in the set of nutrients available during the computation of the flux distribution. accordingly, this approach should provide a solution in which fluxes are computed in such a way that a minimal metabolic adjustment is computed when the model is presented to a different set of nutrients. for sake of clarity, we named this alternative use of moma as nutritional-moma. it is to be noticed that, as data from a t point were not available, our nutritional-moma computations started from phase p. as a result, plots for this set of analyses embed nine time points rather than ten as in fba simulations. regulons and conserved motifs identification the regprecise database was used to inspect the presence and the structure of putative regulons in phtac specific upstream motifs finding was performed using the meme suite. genome-wide regulatory motifs searches were performed using the dminda server on the set of phtac upstream regions present in the door database.  modelling procedures we have recently reconstructed a genome-scale metabolic model of phtac using it for inferring the metabolic adjustments of this bacterium induced by changes in gene expression following a temperature downgrade in this bacterium. here we used this model to investigate the metabolic rearrangements occurring during growth in a nutritionally rich medium. a detailed analysis of the data reported by wilmes et al. allowed the identification of twelve distinct phases in the growth of phtac on peptone medium, each of them corresponding to a time step of one hour. as during the last couple of phases almost no nutrients uptake was recorded, we limited our analyses to the first ten phases. for each of these phases we identified the specific uptake rates of the different compounds present in the growth medium and/or switches in the use of the available c-sources. these values were used as input for ten different fba simulations to derive the most likely fluxes distribution in the phtac metabolic model in each time step. this allowed a system-level characterization of the metabolic changes occurring after variations in the usage of the different carbon sources. a schematic representation of the nutrients provided to the model in each of the time steps is reported in fig. 1a.fig. summary of phtac genome-scale reprogramming following nutrients switching. a. the nutrients provided to the model in each different growth phase according to b. heat map with log values of fluxes across all the phases. c. number of flux carrying reactions in each growth phase. d. number of flux-changing reactions in each growth phase. the dashed line represents the average number of reactions carrying flux over all time points. e. number of reactions whose flux is predicted to increase and decrease following each shift in the nutrients provided; the black line accounts for those reactions whose flux is predicted to decrease not as an effect of an imposed reduced growth rate during simulations. f. functional annotation of reactions varying their flux across all the phases prediction of core and switching reactions first, we evaluated the predicted flux for each reaction in the model, across all the time points. overall, we found that reactions did not carry flux in our model in any of the growth phases; these reactions are probably not essential to sustain phtac125’s growth in the simulated conditions and were discarded from the following analyses. conversely, reactions were predicted to carry flux in at least one time point and were considered for further analyses. on average, excluding the last time point in which all the reactions are predicted to be turned off, reactions are predicted to carry flux for each of the analysed phases. interestingly, this value is quite constant throughout the time points, revealing that, according to our simulations, the number of reactions necessary to support phtac growth is predicted to be somehow independent from the carbon sources used as c and energy sources. despite this apparently constant trend in the number of metabolic reactions employed by phtac in a complex medium, the shift among the different growth phases is presumably characterized by a relatively high number of reactions that are predicted to change their flux. on average, reactions display at least a single variation in their predicted flux across all the time points. however, in this case, a greater variability across all the time points is observed in our simulations. gene-protein-reaction rules of our model indicate that this latter set of reactions is encoded, on average, by genes. interestingly this number resembles that found when modelling the metabolic switch of s. coelicolor and represents around 12% of the phtac coding sequences. according to our simulations, switching among the available carbon sources in a complex medium may have an impact on the overall metabolic network of phtac with more than half of flux-carrying reactions influenced by a change in the utilized c-source or in its uptake rate. more in detail, we found that, on average, reactions are predicted to decrease their flux throughout the growth period, whereas are predicted to increase it. however, the number of reactions whose flux decreases might be biased due to a general decrease, across the time steps, of the growth rate. more specifically, this systemic bias is related to the amino acid uptake rates, which have been derived from the physiological data by wilmes et al., whose constant decrease lead, for some reactions, to a flux reduction in each time step. thus, we adjusted the set of “decreasing reactions”, by removing those for which we observed a consistent decrease for all the growth phases. although the normalization did not affect the general trend, in that the normalized set and the not normalized one have similar trends, the number of decreasing reactions is comparable to that of the increasing reactions. furthermore, we also computed a normalized flux distribution for each of the modelled growth phases, expressing them as a fraction of the predicted growth rate and evaluated whether this procedure led to different results. results revealed no major differences in the overall trends compared to the original calculation of fluxes distribution, suggesting that our results hold true regardless of the normalization procedure adopted to account for the different growth rates predicted by the model across the growth period. not all phtac metabolic pathways are impacted by this switching of nutrients in the same manner, according to our simulations. figure 2a shows the hypothetical number of flux-carrying reactions for five pathways, i.e. tca cycle, lys biosynthesis, glu metabolism, val, leu and ile biosynthesis and degradation. tca cycle, for example, displays and increasing trend in the number of flux-carrying reactions according to our modelling framework; this is consistent with the exhaustion of amino acids whose degradation provides important tca cycle intermediates, i.e. oxaloacetate, fumarate and 2-oxo-glutarate and, consequently, with the necessity to activate those reactions leading to the biosynthesis of such compounds. conversely, lys and glu metabolic routes display an overall constant trend, with a similar number of active reactions across the different simulated growth phases. this is in line with i) the necessity to use of the lysine biosynthetic route to synthesize diaminopimelic acid, an essential component of bacterial cell wall and ii) with the importance of glu metabolism for phtac. finally, val, leu and ile biosynthesis and degradation pathways display an opposite trend one another. intuitively, this reflects the necessity to synthesize these molecules in the first part of the growth phase and the necessity to catabolise them once phtac is using those amino acids as carbon sources, respectively.fig. changes in the central metabolism of phtac a. the number of active reactions for five major pathways across all the time points is shown. b. a simplified representation of the interconnections in the central metabolism of phtac dashed lines indicate the presence of more than one reaction between the connected compounds. modified from we next evaluated whether changes in the kind and number of utilized substrates across all the time points was also parallel to major structural changes at the whole metabolic network level. we used flux variability analysis to perform a comprehensive exploration of alternate optimal routes in phtac metabolic network, in each of the predicted growth phases. this analysis revealed that the number of “flexible reactions” is predicted to remain constant throughout all the growth phases, suggesting that the plasticity of phtac metabolic network is not particularly dependent on the set of metabolized substrates. this suggests the presence of an efficient adaptation of phtac125’s metabolic network to face a nutritionally rich but probably variable environment. overall, according to our model, the central phases of the growth curve are characterized by the largest set of reactions changing their flux, consistently with the major changes in the set of utilized nutrients observed experimentally by wilmes et al.. the first important switch in the panel of utilized compounds occurs between phase and. according to experimental evidence from wilmes et al., in this time point phtac stops using ser, asn and asp as c and n sources, relying only on glutamate for sustaining growth. the importance of glutamate was previously observed from fed-batch cultivation experiments as it was the most strongly metabolised amino acid in all growth experiments. also, glu is the amino acid allowing the fastest growth rate among all the amino acids when growth in twenty minimal growth media  * h) was simulated. during p phtac relies only on glutamate for sustaining growth and needs to utilize this compound to derive its main building blocks. l-glutamate:nadp+ oxidoreductas is predicted to be a key reaction in this phase as it allows the conversion of glutamate to 2-oxoglutarate to be used as energy source feeding the tca cycle. more in general, our modelling framework predicts a deep re-programming at the whole metabolic network level following this shift, with reactions changing their predicted flux; more in detail, reactions are predicted to decrease their flux, whereas are predicted to carry an increased flux. the remaining reactions are predicted to operate in the reverse direction compared to the previous growth phase following the switch. among the first set of reactions, are predicted to be turned on in respect to the previous phase, whereas are predicted to be turned completely off. this latter set of reactions includes, as it might be expected, the transporters of serine, asparagine and aspartate, i.e. those compounds that are no longer available in the medium and that cannot be internalized anymore. the cell faces the absence of some of these compounds by redirecting of the glutamate available in the medium to their synthesis. indeed, for example, the reaction encoded by l-aspartate:l-glutamine amido-ligase leading to the formation of asparagine from gln is predicted to double its flux after this metabolic transition. the exhaustion of asp is probably compensated by reversing reaction r, generating asp from glu. furthermore, our simulation indicated a drop to no flux in the reaction encoded by threonine dehydratase and catalysing the conversion of thr to production of 2-oxobutanoate, a precursor of val, leu and ile. a decrease in the availability of asp might impair the biosynthesis of 2-oxobutanoate following this pathway. our model predicts that phtac faces this perturbation by reversing the flux of reaction r, leading to the production of 2-oxobutanoate from o-succinyl-l-homoserine. after the following switch, wilmes et al. observed that phtac starts using ala, leu, gly and tyr. according to our modelling framework, this involves a predicted change in the flux of reactions, with and of them decreasing, increasing or changing the direction of their flux, respectively. consistently with this new set of amino acids, the model predicts that most of the biosynthetic routes leading to the production of such compounds result turned “off”. this is the case of the path leading to the synthesis of leu in the branched-chain amino acid biosynthetic route. the same occurs for ala biosynthesis and l-aspartate 4-carboxy-lyase are predicted to be turned off) and lys biosyntheses. interestingly, in the latter case, only the final step of the whole pathway is predicted to carry no flux, consistently with the observation that the rest of the pathway is required for the synthesis of meso- -diaminopimelate, an essential precursor for peptidoglycan assembly. amino acid biosynthetic routes are not the only pathways affected by this metabolic switch, according to our simulated growth model. two tca cycle reactions are predicted to be turned on following the utilization of ala, leu, gly and tyr, i.e. those leading to the production of 2-oxoglutarate from isocitrate. this raises the intriguing question on the source of 2-oxoglutarate in the previous growth phases given that those tca reactions were predicted to carry no flux. a possible explanation is provided by the observation that 2-oxoglutarate can be obtained from the carbon skeletons of several five-carbon amino acids through a first conversion into glu, which is then oxidatively deaminated by glutamate dehydrogenase to yield α-ketoglutarate. however, the presence of glu in the medium would not favour this solution as the major fraction of glu can be obtained without additional energy expenses. alternatively, 2-oxoglutarate can be obtained from the conversion of glu to tyr. in our model this reaction is predicted to be “on” during p thus allowing the synthesis of 2-oxoglutarate. however, starting to use tyr present in the medium during p might cause reaction r to be turned off and, consequently, the necessity to synthesize 2-oxoglutarate from succinate to maintain the functioning of tca cycle. similarly, the production of fumarate from malate is predicted to carry no flux in the shift to p. again, this is consistent with the degradation of amino acids as a major source of important metabolic intermediates, as fumarate can be obtained from the degradation of aromatic amino acids as tyr, available in the medium during p and then tunnelled into the tca cycle. the time-resolved growth data from wilmes et al. show that, after other minor transition in which few additional amino acids are degraded, the final switch involved the utilization of histidine as the only c source and causes a major drop in the growth rate of phtac. in our model, this corresponds to the highest number of flux-changing reactions and thus, to the deepest reprogramming encountered by phtac during this growth curve. as it might be expected, histidine biosynthetic reactions are predicted to carry zero flux following this transition. despite the entire pathway is predicted to carry no flux, aicar -5-amino-4-imidazolecarboxamide), one of the intermediates of histidine biosynthesis and a crucial precursor in purine metabolism, might still be synthesized through reaction r. this compound, thanks to the flux predicted to be carried in this phase by reaction r leads to the synthesis of 5-phospho-alpha-d-ribose 1-diphosphate from d-ribose 5-phosphate. furthermore, the switch to the utilization of his as the sole c-source is predicted to be responsible for the stop of the catabolic routes of the previously degraded amino acids and, consequently, in the stop of the production of important cellular intermediates. according to our simulation, this causes the re-activation of key metabolic reactions that cannot rely on many degradation intermediates as in the previous phase. these include those belonging to the tca cycle and purine metabolism. finally, the reactions involved in the conversion of his into glu are predicted to be turned on following this final transition, allowing the production of glutamate and, from these, all the major pathways found to carry flux also in phase p. a sub-optimal objective function predicts alternative fluxes distribution it is worth noticing that all the simulations described up to now were conducted under the fba canonical assumption of biomass optimality, i.e. assuming that all metabolic fluxes in the cell are geared towards the production of biomass in each moment of the growth curve. however, situations in which cells invest substantial resources in alternative metabolic processes at the expense of optimal growth rates have been analysed and described quite extensively. according to this scenario, cells may invest substantial resources in a specific metabolic process at the expense of optimal growth, this being reflected by sub-optimal flux distributions. this situation might be observed, for example, when microbes allocate energetic resources for anticipating changing environmental conditions at the expense of optimal growth. also, when cells are exposed to nutrients fluctuations they might respond with a minimal metabolic adjustment, to avoid the waste of protein synthesis and degradation necessary to reprogram the entire metabolic network and to simultaneously achieve two objectives, i.e. rapid and minimal adjustments. this latter scenario may indeed resemble the actual competition for nutrients that emerges in natural environments and that has been proposed to be reflected by the sequential uptake of nutrients. to explore this alternative scenario, we have accounted for the possibility that the phtac biomass objective function could not be fully optimized but, instead, in a near-optimal or sub-optimal state. since fba classically assumes the optimization of the biomass production flux to compute the most likely fluxes distribution inside the cell, it cannot per se provide hints concerning alternative fluxes distribution. for this reason, we have used moma to study the hypothetical fluxes distribution when minimizing the metabolic adjustments required at each transition of the entire growth period analysed. this differs from the canonical formulation of moma in which the effects of a gene knock-out are evaluated by providing an approximate solution for a sub-optimal growth flux state, nearest in flux distribution to the unperturbed state. applying this modelling strategy to our study case revealed that the choice of the optimization criterion has a great influence on the predicted fluxes distribution, both in terms of the fraction of the entire network required to sustain phtac growth and on the set of active pathways. first, a higher number of active reactions compared to the original fba predictions were predicted for each of the growth phases, when the implemented nutritional-moma approach was used. this gap is even more evident in the first seven growth phases, whereas the difference between the two approaches becomes negligible in the last three phases. however, despite the two approaches predict a similar number of active reactions in these two final time points, the two sets of reactions appear to be quite different.fig. comparison between nutritional-moma and fba. here we show the number of predicted flux carrying reactions in each growth phase for fba and nutritional-moma optimization on the phtac model. also, the number of shared reactions identified by the two approaches is shown the fact that nutritional-moma always predicts a higher number of active reactions might be explained by considering that to optimize the model using fba identifies the flux distribution maximizing the objective function, regardless of the metabolic states in the previous time steps. conversely, nutritional-moma will seek for the solution that is closer to the one of each previous time step. therefore, the results of moma will likely include a higher number of reactions in the model to produce biomass, since it will minimize the changes with respect to the previous time steps, while activating novel reactions to cope with the changes in the nutrient composition. we then investigated whether these differences in the sets of predicted flux-carrying reactions between the two approaches were involved in specific pathways or, rather, spanned over a larger part of the entire phtac metabolic network. thus, for each metabolic pathway in the model, we computed the fraction of reactions predicted to be active by nutritional-moma, fba and both methods. results of this analysis revealed that the choice of the optimization method impacts the predicted fluxes distribution not only for what concerns the activity of the peripheral pathways, i.e. those pathways that start the degradation of amino acids and then tunnel them into the central metabolism. indeed, we found many central processes in which the proportion of active reactions is specific for each of the two optimization strategy. notably, for the first six growth phases, most of the reactions predicted to be active under fba optimizations are carrying flux also adopting the nutritional-moma approach. in some cases, however, entire pathways are differently predicted to be active/inactive by the two approaches. the last two phases of the growth nutritional-moma and fba predict a similar number of active reactions but, from a functional viewpoint, deep differences exist in that entire pathways are predicted to be active only under a specific optimization method.fig. functional differences between nutritional-moma and fba predictions. here we show the proportion of reactions predicted to be active by nutritional-moma, fba and both methods for each main functional category represented in the phtac reconstruction to summarize, using our moma-based approach we identified a higher number of reactions predicted to be carrying flux in respect to the fba-based optimization. the nutritional-moma approach predicted larger active metabolic networks throughout the analysed phases possessing, on average, about flux-carrying reactions more than the one simulated with fba. furthermore, despite most of the active reactions identified through fba were also identified by the nutritional-moma approach, in some cases these predictions differed significantly from a functional viewpoint. at present, further experimental evidences are needed to shed light on the real number and function of reactions used by phtac in each time point and, in other words, to infer how far from the actual fluxes distribution our in silico predictions are.fig. flux correlation analysis. heatmap accounting for the pearson correlation of all the flux difference vectors across all the time points. the metabolic process of each reaction is also reported flux correlation analysis identifies functionally associated genes in this section we analyse more in depth the flux matrix. this matrix contains the fba-predicted flux of the reactions for each of the time steps. as such, it allows capturing the co-variation of reactions in the model; in other words, we identified as paired reactions those whose predicted fluxes displayed a similar trend throughout the time points. these, in turn, may represent functional partners in the cell and/or may belong to the same metabolic pathway/module. we represented co-varying reactions in the form of a heatmap as shown in fig.  in most cases, the small clusters of reactions showing high pearson correlation values are involved in the same metabolic process. this is compatible with a scenario in which, following a nutrient switch, entire metabolic pathways are activated to face the novel environmental condition. nevertheless, grouping reactions with a lower but still significant threshold, clusters start to embed reactions belonging to other metabolic processes. these, in turn, may represent previously undetected functional associations between genes and/or entire pathways. overall, our method led to the identified different clusters fig. 6a, comprising reactions. according to the gpr of our model, these reactions were encoded by genes.fig. a. co-varying reactions clusters identification. common flux trends for the reactions embedded in each of the clusters. b. the distribution of string combined scores among all the genes embedded in each cluster of genes and the number of genes embedded by each cluster. the grey line represents the median of the combined score computed for each possible pair of genes in the model we analysed these clusters to answer two biological questions, namely whether the embedded genes i) are de facto functionally related and ii) share conserved upstream motifs. groups of co-varying genes display higher combined evidence scores to address the first point, we exploited the combined evidence score provided by the string database. as shown in fig. 6b, the median of the combined score for the genes embedded in the same cluster is, in most cases, well above the median of all the possible combinations of the genes embedded in the model, suggesting that flux coupling in the model occurs among bona fide functional patterns. accordingly, this result corroborates the capability of our modelling framework to identify functional metabolic modules and the functional associations of their encoding genes in phtac interestingly, in some cases apparently un-related genes seem to be embedded in the clusters. despite these instances may represent erroneous predictions of our model, they might also suggest the presence of still undetected functional associations and/or the use of common pools of chemical intermediates by the corresponding reactions. we here present a description of the functional relationships retrieved for genes belonging embedded in three of these clusters. a detailed list of the reactions and genes embedded in all the clusters is provided as additional file to allow further experimentation and/or in silico analyses.fig. sample string clusters. evidence network, co-expression instances and co-occurrence patterns for clusters, and. red asterisks in g indicate those genes known to belong to the same regulon according to regprecise database cluster this cluster is composed of seven reactions encoded by genes. three of these reactions are predicted to be catalysed by enzymes encoded by multiple genes, i.e. cytochrome-c oxidase encoded by pshaa2869- ferrocytochrome-coxygen oxidoreductase encoded by pshaa1842- and atp synthase encoded by pshaa3007- according to the gpr rules of our model, the simultaneous presence of all the corresponding genes is required to carry flux of each of these reactions. the functional association network obtained probing the string database with all the genes of cluster is reported in fig. 7a. genes encoding proteins involved in the same enzymatic reaction appear to be highly interconnected among themselves with evidences that included co-expression, gene fusions, and experimental evidences. this is not surprising, as it might be expected that genes encoding proteins that participate in the same enzymatic complex should display coupled regulation and co-occurrence patterns in other organisms. links, however, are also found among members of different enzymatic complexes, consistent with the fact that these reactions belong to the overall process of oxidative phosphorylation. the ndk gene is also found connected to one of the sub-cluster fig. 7a). this protein plays a major role in the synthesis of nucleoside triphosphates and is involved in purine and pyrimidine metabolism. the coupling of this reaction’s flux with those belonging to oxidative phosphorylation can thus be explained by their common final products, as optimizing the model for biomass production would probably force the flux through both these pathways. finally, the fact that orthologs of these genes are found to be co-expressed and located in physical proximity in other genomes support the functional association predicted by our model. cluster we found eight reactions embedded in this cluster, involved in the following metabolic pathways: glycerolipid metabolism, pentose phosphate pathway and glycolysis. half of the reactions of this cluster are involved in the conversion of glycerol into glycerate 3-pshosphate that, in turn, is connected to glycolysis/gluconeogenesis. as shown by the string analysis output, many evidences exist of a functional association among the genes of this cluster, in most cases including instances of co-expression and physical proximity fig. 7e and f, respectively). two genes of the cluster remained disconnected in the string network, both encoding a glycerol kinase and assigned to the same reaction in our metabolic reconstruction. accordingly, despite no previously detected functional associations were retrieved between these two genes and the other embedded in the same cluster, their presence could be accounted for by the use of common precursors with the reactions encoded by cluster genes. indeed, fluctuations in the availability of such compound during the simulated growth transitions explains the coupling of such reactions. cluster this cluster embeds nine reactions encoded by genes involved in arg biosynthetic pathway argah, b, c, d, e and in the metabolism of pro ). the string-based evidence network in fig. 7g revealed hints on the functional associations existing among the genes of the cluster including co-occurrence of orthologs of these genes in other genomes, as well as co-expression and conserved genomic neighbourhoods fig. 7h and i, respectively). the genes embedded in this cluster encode enzymes involved in the the formation of arg and pro from glu, using asp and then orn as intermediates. these reactions appear to carry flux especially in the first three phases of the growth, that is when aspartate is present in the medium. afterwards, a decrease in the flux carried by these reactions is observed, followed by an overall constant flux trend in the next transitions. accordingly, the correlation among the fluxes of these reactions could be explained by the need to turn these pathways “on” when glutamate and required intermediates are present in the growth medium, respectively. a gene encoding an acetyl-coa synthetase is also embedded in this cluster, despite it has not connections with arginine and proline metabolism-related genes. however, its presence might reflect the necessity to synthesize acetyl-coa from acetate during the first four growth phases identified, i.e. in the absence of ketogenic amino acids. once phtac starts metabolising the second set of amino acids, acetyl-coa can be obtained from the degradation of specific amino acids and thus the reaction encoded by acsa is predicted to stop carrying flux. analysis of upstream sequences identifies putative regulation mechanisms in the previous section we reported a list of gene sets which are, presumably, functional partners and whose activity is concerted during the different metabolic switches. if this holds, it is reasonable to hypothesize that some of these clusters might embed co-regulated genes that, in turn, may share conserved upstream motifs. thus, to address this point we applied an ad hoc computational pipeline to analyse the upstream regions of the genes embedded in each cluster. this resulted in the identification of motifs that are shared in the upstream regions of all the genes for each cluster. upstream regions were identified using door. these have been annotated using three different motif databases, to identify motifs with putatively related functions and that might be recognized by putatively related transcription factors. the whole set of annotations has been manually curated in a conservative way to obtain high-confidence motif annotations. finally, the obtained putative tfs have been cross-validated by comparing their associated functions with those related to the genes embedded in the clusters. the whole set of data related to the identified motifs and their annotations, for each database, is provided as additional file on average, we identified for each cluster a large number of annotations according to the different databases. however, after manual curation, the number of annotations was reduced to divided in clusters. most of the remaining annotations come from prodoric, while the other databases have a minor contribution of high-quality annotations. for each cluster, we eventually compared the putative tfs with the functions of the genes embedded in the cluster, using both literature information and the regprecise database. with few exceptions that will be described in details, we found no correlation between the putative tfs and the regulation/biological function associated to the genes embedded in the clusters, suggesting that the majority of the co-varying groups share a common regulation but are regulated by a number of different tfs. we are describing below the clusters for which we found agreement between the predicted tfs and the function of the embedded genes.table main features of the putatively co-regulated clusters found during flux-correlation analysis. in this table, for each cluster, we report its number, the name of the regulator identified, the genes embedded in it and the conserved motif found upstream of its genes cluster this cluster, as described previously, is involved in the metabolism of arg and pro. a search on the regprecise regulon database revealed that, at least in the case of arg, the genes embedded in the cluster of fig. 7g are part of a common regulon. the only exception is represented by argd which, according to the regprecise database is not part of the arginine regulon in phtac and the other genes involved in the biosynthesis of pro. however, both the computational approach described above and a manual search for putative transcription factor binding sites upstream of these genes detected the already known conserved motif shared by the other arg genes of the regulon. this suggests that the all the genes are part of the arginine regulon in phtac more in general, this finding confirms that, starting from our modelling outcomes it is possible to infer biologically consistent patterns of gene co-expression. cluster this cluster includes genes from the sdh operon, encoding the succinate dehydrogenases system, as well as other enzymes involved in energy conversion, such as other dehydrogenases and other enzymes with a role in nucleotide modifications. our computational pipeline predicted two different regulator genes associated with this cluster, i.e. ccpa and psra. these regulators encode proteins are involved in the carbon and fatty acid metabolism, respectively. comparing the functional role of the tfs with the cluster annotation, we can easily relate the set of genes encoding dehydrogenases with ccpa and psra, which is also confirmed by regprecise database. the other three genes are indirectly involved in the oxidative metabolism, in that they are involved in the metabolism of fad. cluster this cluster includes five genes, with three of them involved in the synthesis of folate, while the remaining two encode a galactokinase and a galactose-1-epimerase. all these genes are involved in the metabolism of nucleotides, since folate is necessary to synthesize tdp sugars from udp, while the other genes encode enzymes that are involved in the synthesis of udp. the analysis of the upstream sequence of these genes revealed a conserved motif similar to the galr binding site, thus suggesting that these genes are being co-regulated.  here we have used constraint-based metabolic modelling to provide a time-resolved, systems-level scheme of phtac metabolic re-programming following nutrients switching in a nutritionally complex medium. such features have been analysed using the metabolic model and growth data of the antarctic bacteria p. haloplanktis tac previous experimental tests revealed a number of nutrients switches in this microorganism when grown in a complex medium, consistent with the fact that this antarctic organism is adapted to fast growth in a fairly rich environment. indeed, sequential uptake of nutrients is thought to emerge when competition for nutrients is present. our modelling framework identified the central phase of the growth curve as a probable key reprogramming point, with more than reactions predicted to adjust their flux in these time points. this corresponds to the exhaustion of most of the first metabolized c sources, a time step in which only glu is used as the sole c source, and the final part of the growth curve in which phtac stops using glu and relies on a completely different set of nutrients. our model highlights the occurrence of such an adaption and the need for reprogramming a large set of reactions to maintain an efficient metabolic network. a similar scenario is observed at the end of the growth, when his is used as the sole c source. this transition is the most demanding, requiring a change in the predicted flux of almost reactions. taken together, our simulation indicates the presence of an almost constant number of reactions that are required to sustain life across all the time points examined. interestingly, this set of reactions resembles, in percentage, the size of the minimal metabolic network predicted to be active in e. coli. to maintain such a constant trend, however, a deep reprogramming of the whole cellular metabolism appears to be necessary during the entire growth period, as up to more than biochemical reaction, display at least one change in the predicted carried flux. according to our simulations, these changes do not only involve peripheral metabolic pathways but also a number of central pathways, e.g. tca cycle, glycolysis and pp pathway. tca cycle, for example, displays an increasing trend in the number of flux-carrying reactions parallel to the exhaustion of the amino acids whose degradation provide key tca cycle intermediates. this is in line with the assumption that fluxes distribution within the cell are influenced by the entry-point of the a given c-source into the metabolic network. indeed, the c sources provided in our simulations have different distances from the phtac central metabolism and thus are expected to cause a re-wiring of an important fraction of the network. this, in turn, suggests the presence of an efficient metabolic reprogramming machinery to continuously and promptly adapt to this nutritionally changing environment and/or to the exhaustion of the preferred carbon source. modelling and dividing the growth curve in discrete time points has allowed us to infer common trends in the predicted flux patterns of the reactions of the model. we have shown that such coupled patterns likely correspond to reactions and enzymes that work in a concerted fashion and that, in some cases, represent functional partners and/or be encoded by co-regulated genes. the modelling framework we have set up here has allowed gaining insights on the regulation of phtac metabolic genes and, as in the case of the arg regulon, to expand the current knowledge on commonly regulated genes. furthermore, we have assembled a dataset of putatively co-regulated genes that can be used for further manual curation and/or subject to experimental validation. despite fba-based predictions led to the identification of biologically consistent trends, the nutritional-moma approach we have implemented here suggests that care needs to be taken when choosing the objective function for constraint-based metabolic modelling. indeed, when we accounted for the possibility that cells could be in a near-optimal or sub-optimal state the two modelling frameworks predicted quite a different topological and functional reprogramming of phtac cells. as it is currently not possible to unambiguously discern among these two alternative solutions provided by fba and nutritional-moma, further experimental evidence of phtac cells grown in complex medium will allow deriving a clearer picture of those pathways that are really active during the growth and, consequently, which of the two approaches outperforms the other. additional files additional file 1: the metabolic reconstruction and the matlab codes used in this work. additional file 2: additional analyses on the phtac metabolic reconstruction showing i) the effect of expressing fluxes as a fraction of the biomass production flux and ii) number of active reactions in each phtac metabolic pathway following the nutrients switching. additional file 3: the list of the reactions and genes embedded in all the clusters identified through flux correlation analysis. additional file 4: the whole set of data related to the identified motifs and their annotations, for each probed database. abbreviations fbaflux balance analysis fvaflux variability analysis momaminimization of metabolic adjustment phtac125 pseudoalteromonas haloplanktis tac125 tcatricarboxylic acid tftranscription factor 
in their natural environment, bacteria face a wide range of environmental conditions that change over time and that impose continuous rearrangements at all the cellular levels. we used multi-step constraint-based metabolic modelling to simulate the growth in a complex medium over several time steps of the antarctic model organism pseudoalteromonas haloplanktis tac as each of these phases is characterized by a specific set of amino acids to be used as carbon and energy source our modelling framework describes the major consequences of nutrients switching at the system level. here we provide a time-resolved, systems-level scheme of phtac metabolic re-wiring as a consequence of carbon source switching in a nutritionally complex medium. our analyses suggest the presence of a potential efficient metabolic reprogramming machinery to continuously and promptly adapt to this nutritionally changing environment, consistent with adaptation to fast growth in a fairly, but probably inconstant and highly competitive, environment. the model predicts that a deep metabolic reprogramming might be required to achieve optimal biomass production in different stages of growth, with at least half of the cellular metabolic network involved. growth simulations under this scenario revealed the deep impact of choosing among alternative objective functions on the resulting predictions of fluxes distribution. a systemic re-organization of the overall microbial metabolic network in response to a variation in the composition/concentration of the surrounding nutrients has been suggested, although the range and the entity of such modifications in organisms other than a few model microbes has been scarcely described up to now. when facing a nutritionally rich environment, for example, microbes first use the preferred compound and only later start metabolizing the other one. additionally, we show that our modelling framework is able to capture metabolic functional association and/or common regulatory features of the genes embedded in our reconstruction.
 in their natural environment, bacteria face a wide range of environmental conditions that change over time and that impose continuous rearrangements at all the cellular levels. when facing a nutritionally rich environment, for example, microbes first use the preferred compound and only later start metabolizing the other one. a systemic re-organization of the overall microbial metabolic network in response to a variation in the composition/concentration of the surrounding nutrients has been suggested, although the range and the entity of such modifications in organisms other than a few model microbes has been scarcely described up to now.  we used multi-step constraint-based metabolic modelling to simulate the growth in a complex medium over several time steps of the antarctic model organism pseudoalteromonas haloplanktis tac as each of these phases is characterized by a specific set of amino acids to be used as carbon and energy source our modelling framework describes the major consequences of nutrients switching at the system level. the model predicts that a deep metabolic reprogramming might be required to achieve optimal biomass production in different stages of growth, with at least half of the cellular metabolic network involved. additionally, we show that our modelling framework is able to capture metabolic functional association and/or common regulatory features of the genes embedded in our reconstruction. finally, to explore the possibility of a sub-optimal biomass objective function we have implemented a moma-based approach and compared the outcomes with those obtained with flux balance analysis. growth simulations under this scenario revealed the deep impact of choosing among alternative objective functions on the resulting predictions of fluxes distribution.  here we provide a time-resolved, systems-level scheme of phtac metabolic re-wiring as a consequence of carbon source switching in a nutritionally complex medium. our analyses suggest the presence of a potential efficient metabolic reprogramming machinery to continuously and promptly adapt to this nutritionally changing environment, consistent with adaptation to fast growth in a fairly, but probably inconstant and highly competitive, environment. also, we show i) how functional partnership and co-regulation features can be predicted by integrating multi-step constraint-based metabolic modelling with fed-batch growth data and ii) that performing simulations under a sub-optimal objective function may lead to different flux distributions in respect to canonical fba. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords flux balance analysispseudoalteromonas haloplanktis tac125antarctic bacteriametabolic modellingpnrapnra13-b pnra13-az fondi marco fani renato issue-copyright-statement© the author 2016 