Rational system modeling and comprehensive system analysis can serve as prior guidelines for understanding and deducing biological mechanisms. We can retrieve quantitative knowledge for assessing an organism’s metabolic capacity and use this knowledge for in-lab experiments to develop new strains with advantageous productivity [3, 4], or optimizing the cultivation process of existing strains [2].

### Model improvements

Since many studies related to *C. acetobutylicum* ABE pathway have been reported (including parameter values and the rate equation formulas [14, 21]), kinetic modeling of the ABE pathway becomes feasible and enables us to simulate the system dynamics. Nevertheless, the previous kinetic model of *C. acetobutylicum* ABE process (Shinto *et al*, 2007) has several drawbacks as described in earlier context. To overcome these drawbacks, we have established a new model featured with three improvements over the previous one.

First, we have incorporated the key metabolite BuP, reflecting the relevant biological events that are specific to ABE kinetics [15–18]. The correspondence between BuP concentration climax and solventogenesis onset is not merely a natural consequence of the fact that BuP is the intermediate between butyryl-CoA and butyrate. There are implications on the genetic level as stated in Zhao *et al*’s study [16–19]. There are many important solventogenic genes, such as *adh*E1 (CAP0162), *adh*E2 (CAP0035 ), *ctf*A (CAP0163), *ctf*B (CAP0164), *adc* (CAP0165), *bdh*A (CAC3298), *bdh*B (CAC3299), etc., having expression profiles that show a strictly correlated pattern with the kinetics of BuP. Although the detailed mechanism of how BuP acts to regulate ABE process has not been very clear yet, its functional importance has been experimentally confirmed [16–19]. Our new model accounts for this knowledge and is successful in representing the phenomenon.

Second, we describe the regulatory effects of complex factors using a time division pattern. In Shinto’s model, the metabolic regulation beyond the level of substrate/product inhibition/activation is simply defined as the input of glucose. The shut-downs of several acidogenic/solventogenic enzymes (like PTB, BDH, etc) are solely due to the insufficiency of glucose. However, various evidences indicate that even with sufficient supply of glucose, the acidogenic enzymes are still shut down in the solventogenic phase, and the solventogenic enzymes are necessarily inactivated at the beginning of the acidogenic phase [15, 22, 23]. Therefore, the metabolic regulations are not of the simple pattern as Shinto suggested, but a significant 2-phase mode is shown (acids are generated during the earlier phase and solvents are generated during the latter one). In our work, this mode is approximated by considering endogenous enzyme activity variations, assuming enzymes are regulated by many factors (e.g. transcription control) to exhibit different activity levels to fulfill conditional system requirements of different periods. This assumption is equivalent to extending the application of biochemical system theory (BST). In BST, which is based on *in vitro* experiments, enzyme concentrations and endogenous enzyme activity levels are constant by default. Hence kinetic models based on BST are rigorously suitable for chemical simulations but may not entirely appropriate for *in vivo* conditions. Under *in vivo* conditions, the rate of a reaction does not solely depend on substrate/product concentrations, because the endogenous enzyme activity itself is regulated by many factors and its variation in turn affects the reaction rate [15–17]. Our model divides time into a set of periods according to the enzymes’ activity variations, allowing enzyme activities to vary throughout these periods.

Third, we introduce the “enzyme activity coefficient (EAC)” to quantify endogenous enzyme activity variations caused by metabolic regulations (see section “Methods” for EAC’s definition). For the quantification of enzyme activity curves, numerical interpolation (e.g. Lagrange, Legendre, etc.) should have been employed as to obtain fully continuous functions. But measurements in activity assays are usually not precise. If the errors are large, interpolation may result in huge errors or mistakes, causing the trouble of overfitting and distorting the original curve profile. On the contrary, the computation of EAC leaves the error just as the original error. Hence, using EAC will at least not amplify the error or distort the curve when the measurements are not precise. Moreover, our design of EAC is calculating a ratio instead of the particular value at each time instance, and this allows the error to be divided by a denominator, thus lowering the error level in computation.

### Dynamical simulation and perturbation analysis

After the addition of BuP, 5 unknown parameters are introduced into the system. We have used Genetic Algorithm to estimate their values. In the process of parameter estimation, we used Shinto’s experimental observations of 16 metabolites to formulate the fitness function, but we didn’t employ any information about BuP. And in order to avoid the mistake of reasoning in a circle, we compare our results with observations of another experiment (Zhao *et al.*’s). It turns out that our results are significantly consistent with the observations and have shown some superiority over Shinto’s model in reflecting the kinetics of BuP and butyrate. This indicates that Shinto’s model is well fitted for its own condition but may not be suited well for other conditions. In contrast, our model has more capacity in approaching real cases because of the improvements we have made.

Simulations based on kinetic models can help develop in-lab strategies, thus increasing the success rate of metabolic engineering. In our work, we have simulated thousands of perturbed conditions to detect and assess potential spots that have large influences on butanol production. The magnitude of *in silico* perturbations should not be too large because the system may exhibit alternative activations for other pathways when undergoing substantial fluctuations [24, 25]. When the system is encountering slight perturbations, its overall properties will not change substantially due to biological robustness [25–27]. So it would be fairly assumed that when the perturbation magnitude on enzymatic parameters is 5%, the system will still survive and its functional normality is not interrupted or diverted. In the computation, we have identified an interesting phenomenon that BK’s catalytic capacity exhibits positive influence on butanol production while CoAT has negative influences, as elevating BK activity results in increased *Rd* value and uplifting CoAT activity diminishes the value. And more convincingly, *Rd* decreases when increasing the *V*
_{
max
} values of BK and CoAT at the same time, which means the negative effect of CoAT can balance the positive effect of BK, confirming that CoAT has large effect in impairing butanol production. Based on this discovery, we propose a possible scenario that if more metabolites are received by BK as substrates, the overall acid (butyrate) reassimilation efficiency will be benefited and butanol production is enhanced. And if more metabolites are received by CoAT as substrates, the situation will be on the contrary. It may not seem economical for the bacteria to use the BK-PTB path (see Figure 1) to reassimilate butyrate since running through this path consumes ATP. Nonetheless, based on our computation results and biochemical knowledge, we raise a predictive explanation for the underlying mechanism: in acidogenic phase, the metabolic flux actually runs in the direction of PTB-BK (confirmed by both our computation of kinetic profiles and experimental literature [5, 16, 21]), thus this path generates ATP for the growth of the bacteria; when the bacteria enters solventogenic phase, it doesn’t need to grow and ATP has surplus, these surplus ATPs are utilized to proceed butyrate reassimilation. It’s noteworthy that acids are severely poisonous to bacteria cells and it is a priority for the bacteria to convert acids to other forms (e.g. alcohol). In addition, enhanced butanol production means more acids are converted. Hence, although reactions through BK cost ATP, but so far as BK’s efficiency is concerned, BK is still the preferred enzyme through which the bacteria reassimilates butyrate during solventogenesis. Therefore, the reason why the ATP-costing path BK-PTB is more efficient over the path catalyzed by CoAT (not ATP-costing) in reassimilating butyrate is probably because of responding to severe poison stress, and the energetic basis for this process is the ATP surplus generated during acidogenic phase. Our prediction is equivalent to considering the bacterial cellular behaviour to be related with biological robustness, as supposing that the bacteria is not seeking for its optimality when undergoing stress response, but seeking for sub-optimality. In such case, certain costs or sacrifices are tolerated as long as it can survive (or maintain minimal fluctuation from normality) [24, 26].

In double parameter perturbation tests, we noticed that the net effect of combinatorial perturbation was equal to the sum of effects of individual perturbations, indicating that no crossover or nonlinear amplification originated from perturbations with mild magnitudes. This is probably because when the system is undergoing mild perturbation, it tries to maintain the normal status with minor alterations by means of system robustness. To demonstrate the hypothesis further, we implemented some three-parameter combinatorial perturbation tests. We randomly chose a number of three-parameter triplets and randomly decided their shift directions. For example, if we increased three parameters *V*
_{max14}, *V*
_{max19}, *V*
_{max17} by 5% each and re-computed our model (Equation (1)), we obtained *Rd*=2.14%, which exactly equalled the sum of individual effects of these perturbations. Again, if we increased *K*
_{m15b} and *V*
_{max19} by 5% and decreased *V*
_{max18} by 5%, we obtained *Rd*=2.09%, still equalled to the sum of individual effects. Hence, we raise a hypothetic measure for increased butanol production: By slightly perturbing parameters in suitable directions and with appropriately mild magnitudes, we possibly can obtain a metabolic phenotype that can have amplified butanol production, and the strain can steadily and safely survive as well. The amplification magnitudes of multiple parameter perturbations can be much greater than those in single parameter perturbations, if adequately many parameters are manipulated appropriately. Meanwhile, from an engineering point of view, multiple spot modifications can make the risks of system fluctuations or external impulsions more distributive than in the case that all alterations are concentrated on a single spot. Hence, this strategy provides a way that can make a more stable high-production system. But this strategy requires high-precision genetic manipulation.

### Significance

Traditional kinetic models cannot accommodate complex metabolic regulation effects (e.g. gene transcriptional control). Hence previous integrative modeling approaches for metabolic system are mainly based on the FBA method, in which the gene transcription regulations are described by Boolean logic and the metabolic level is expressed by flux balance equations. Since FBA based methods and Boolean logic cannot adequately reflect system dynamics, we have developed a new model as an attempt towards solving the problem. Actually, our modeling strategy is equivalent to extending the traditional BST, degenerating complex metabolic regulation effects to a form that is compatible with kinetic models. This strategy provides a way for integrating complex factors and knowledge from multiple levels into the framework of kinetic models. Moreover, our approach of describing metabolic regulation effects with a time division pattern and EAC is extendable. For instance, we can relate the enzyme activities to gene transcriptional level, build a formulism between them, and include the effects of other factors such as impulse and stochasticity. Our modeling method can be generalized and extended to the modeling of other bio-processes.

In this post-genomic era, massive information and experimental data have been accumulated. Therefore, it is important to develop methods or tools that are able to make use of existed information/data and capable of organizing, manipulating and interpreting them more comprehensively [28, 29]. Our work just attempts to serve that goal by integrating existed information from multiple aspects and describing them mathematically. Nevertheless, the usage of “net effects of regulatory factors” in our modeling doesn’t seem to build direct links between the genetic level and metabolic level. But if adequately more information about the regulatory factors on the genetic level is revealed, better formulism can be built to link the two levels and further studies on the control of bacteria cellular systems can be conducted.