Search

The calculation of steady-state metabolite concentrations in metabolic reaction network models is the first step in the sensitivity analysis of a metabolic reaction system described by differential equations. However, this calculation becomes very difficult when the number of differential equations is more than 100. In the present study, therefore, we investigated a calculation procedure for obtaining true steady-state metabolite concentrations both efficiently and accurately even in large-scale network models. For convenience, a linear pathway model composed of a simple Michaelis-Menten rate law and two TCA cycle models were used as case studies. The calculation procedure is as follows: first solve the differential equations by a numerical method for solving initial-value problems until the upper several digits of the calculated values stabilize, and then use these values as initial guesses for a root-finding technique. An intensive investigation indicates that the S-system technique, finding roots in logarithmic space and providing a broader convergence region, is superior to the Newton-Raphson technique, and the algorithm using the S-system technique successfully provides true steady-state values with machine accuracy even with 1,500 differential equations. The complex-step method is also shown to contribute to shortening the calculation time and enhancing the accuracy. The program code has been deposited to https://github.com/BioprocessdesignLab/Steadystateconc .

Mathematical modeling of large-scale metabolic networks usually requires smoothing of metabolite time-series data to account for measurement or biological errors. Accordingly, the accuracy of smoothing curves strongly affects the subsequent estimation of model parameters. Here, an efficient parametric method is proposed for smoothing metabolite time-series data, and its performance is evaluated. To simplify parameter estimation, the method uses S-system-type equations with simple power law-type efflux terms. Iterative calculation using this method was found to readily converge, because parameters are estimated stepwise. Importantly, smoothing curves are determined so that metabolite concentrations satisfy mass balances. Furthermore, the slopes of smoothing curves are useful in estimating parameters, because they are probably close to their true behaviors regardless of errors that may be present in the actual data. Finally, calculations for each differential equation were found to converge in much less than one second if initial parameters are set at appropriate (guessed) values.

We have previously developed the software for calculation of dynamic sensitivities, SoftCADS, in which one can calculate dynamic sensitivities with high accuracy by just setting the differential equations for metabolite concentrations. However, SoftCADS did not always provide calculated values with the machine accuracy of a computer, although a Taylor series method was employed to numerically solve the differential equations. This is because numerical derivatives calculated from an approximate formula were directly used in the derivation of the differential equations for sensitivities from those for metabolite concentrations. The present work therefore attempts to further enhance the performance of SoftCADS, including not only the accuracies of the calculated values but also the calculation time. To overcome the problem, the approximate formula is expanded into a Taylor series in time and the first-term value of the series is replaced by the exact coefficient on the second term of the flux function expanded into a Taylor series in an independent or dependent variable. The result reveals that this replacement certainly provides not only numerical derivatives but also dynamic sensitivities with superhigh accuracies comparable to the machine accuracy, regardless of the degree of stiffness of the differential equations. Moreover, a comparison indicates that the improved SoftCADS shortens the calculation time of the dynamic sensitivities without reducing their accuracies, even when the simplest approximate derivative formula is used.

This paper investigates the accuracy of a matrix method proposed by other researchers to calculate time-varying flux control coefficients (dynamic FCCs) from elasticity coefficients by means of summation and connectivity theorems in the framework of metabolic control analysis. A mathematical model for the fed-batch penicillin V fermentation process is used as a case example for discussion. Calculated results reveal that this method produces significant calculation errors because the theorems are essentially valid only in steady state, although it may provide rough time-transient behaviors of FCCs. Strictly, therefore, dynamic FCCs should be directly calculated from the differential equations for metabolite concentrations and sensitivities.

Dynamic sensitivity analysis has become an important tool to successfully characterize all sorts of biological systems. However, when the analysis is carried out on large scale systems, it becomes imperative to employ a highly accurate computational method in order to obtain reliable values. Furthermore, the preliminary laborious mathematical operations required by current software before the computation of dynamic sensitivities makes it inconvenient for a significant number of unacquainted users. To satisfy these needs, the present work investigates a newly developed algorithm consisting of a combination of Taylor series method that can directly execute Taylor expansions for simultaneous non-linear-differential equations and a simple but highly-accurate numerical differentiation method based on finite-difference formulas. Applications to three examples of biochemical systems indicate that the proposed method makes it possible to compute the dynamic sensitivity values with highly-reliable accuracies and also allows to readily compute them by setting up only the differential equations for metabolite concentrations in the computer program. Also, it is found that the Padé approximation introduced in the Taylor series method shortens the computation time greatly because it stabilizes the computation so that it allows us to use larger stepsizes in the numerical integration. Consequently, the calculated results suggest that the proposed computational method, in addition to being user-friendly, makes it possible to perform dynamic sensitivity analysis in large-scale metabolic reaction systems both efficiently and reliably.

A simple method for the calculation of dynamic logarithmic gains, i.e., normalized sensitivities, is applied to four typical chaotic systems (Lorenz model, Rossler’s equations, double scroll attractor, and Langford’s equations) to examine an effect of the chaotic behavior on the dynamic logarithmic gains. As a result, it is found that the dynamic logarithmic gains increase exponentially while the chaotic behavior is observed. It is also shown that the observation of the dynamic logarithmic gains is useful to accurately identify the generation of the chaos.

The analysis of immobilized enzyme reactions leads to a two-point boundary value problem that requires numerical solutions based on root-finding methods. The standard Newton-Raphson (N-R) algorithm is adequate if the root-finding is initiated close to the true solution; in other cases, large numbers of iterations may be needed before the algorithm converges. An alternative is a root-finding method based on the theory of S-systems. This method turns out to be much more robust and efficient than the N-R algorithm. It converges quadratically, with fewer iterations, and in a wider range of relevant initial values. An investigation of a simpler algebraic equation with a similar mathematical structure demonstrates that the rate of convergence in the S-system method is not very sensitive to the distance between the root and initial value, whereas convergence of the N-R algorithm slows down drastically with increasing distance.

A variable-order, variable-step Taylor-series method in Cartesian space is discussed which makes it possible to solve simultaneous first-order differential equations expressed in GMA-system canonical form with a super high-order accuracy that is comparable to the machine accuracy of the computer. In this method, the stepsize for each differential equation is calculated from the formula derived by setting the ratio of the appropriate two terms of the first few Taylor-series terms at unity, and a minimum value of them is then used as the initial value of the stepsize at each step of the integration and is halved if several criterions are satisfied. Also, when each dependent variable takes a value near zero, the integration is carried out with a very small stepsize. The Taylor-series method constructed here is first applied to the Lotka-Volterra equation to confirm the validity of the calculation algorithm and it is found that the calculated values completely agree with those by the conventional method in which the Taylor-series solution in logarithmic space is used. Second, the differential equation with a solution of cos(t) is solved and the relative error is found to instantaneously increase when the dependent variables take a value near zero and becomes equal to zero or almost zero in the other region. Finally, the differential equation for a two-body problem is solved and it is shown that the accuracies of the calculated values are kept at high levels even after integration over 625 cycles. In conclusion, the Taylor-series method constructed in Cartesian space is considered to have high applicability to simultaneous differential equations and provides high accuracies for the calculated values.

A numerical method for solving an axial dispersion model (two-point boundary value problem) with extremely high-order accuracy is presented. In this method, one first recasts fundamental differential equations into S-system (synergistic and saturable system) canonical form and then solves the resulting set of simultaneous first-order differential equations by the shooting method combined with a variable-order, variable-step Taylor series method. As a result, it is found that over wide ranges of systemic parameters (Peclet number, dimensionless kinetic constant, and reaction order), this method promises numerical solutions with the superhigh-order accuracy that is comparable to the machine accuracy of the computer used. The advantage of the numerical method is also discussed.

To establish a highly-reliable numerical method for calculation of the values of apparent kinetic parameters in various types of immobilized enzyme reactions, the performance of the shooting method combined with a variable-order, variable-step Taylor-series method has been discussed in a series of two papers. In the first paper, every formulas are given and the calculation technique is described. In this second paper, these formulas are used to calculate the values of apparent kinetic parameters for the Thiele modulus varying from 1 to 700. As a result, the calculated values are found to have extremely-high accuracies that are comparable with the machine accuracy of the personal computer used, and the execution time is always within 1 s. On the other hand, the accuracies of the values calculated by the orthogonal collocation method are shown to reduce remarkably with increasing magnitude of the Thiele modulus.

To establish a highly reliable numerical method for calculation of the values of apparent kinetic parameters in various types of immobilized enzyme reactions, the performance of the shooting method combined with a variable-order, variable-step Taylor-series method has been discussed in a series of two papers. In this first paper, every formulae that are necessary to execute the numerical calculation are given and the procedure and technique for the calculation are described. As an example, the effectiveness factor is calculated for the Thiele modulus from 1 to 700 and the calculated values are found to have accuracies that are almost equal to the machine accuracy of the personal computer used. In the Taylor-series method, the step sizes are efficiently estimated from the ratio of the appropriate two of the first few terms in each Taylor-series solution to the relevant simultaneous differential equations, which contributes greatly to reduction in both the loss-of-significance error and the execution time.

Metabolic reaction systems can be modeled easily in terms of S-system type equations if their metabolic maps are available. This study therefore proposes a method for estimating parameters in decoupled S-system equations on the basis of the Newton-Raphson method and elucidates the performance of this estimation method. Parameter estimation from the time-course data of metabolite concentrations reveals that the parameters estimated are highly accurate, indicating that the estimation algorithm has been constructed correctly. The number of iterations is small and the calculation converges in a very short time (usually less than 1s). The method is also applied to time course data with noise and found to estimate parameters efficiently. Results indicate that the present method has the potential to be extended to a method for estimating parameters in large-scale metabolic reaction systems.

Recent development of high-throughput analytical techniques has made it possible to qualitatively identify a number of metabolites simultaneously. Correlation and multivariate analyses such as principal component analysis have been widely used to analyse those data and evaluate correlations among the metabolic profiles. However, these analyses cannot simultaneously carry out identification of metabolic reaction networks and prediction of dynamic behaviour of metabolites in the networks. The present study, therefore, proposes a new approach consisting of a combination of statistical technique and mathematical modelling approach to identify and predict a probable metabolic reaction network from time-series data of metabolite concentrations and simultaneously construct its mathematical model. Firstly, regression functions are fitted to experimental data by the locally estimated scatter plot smoothing method. Secondly, the fitted result is analysed by the bivariate Granger causality test to determine which metabolites cause the change in other metabolite concentrations and remove less related metabolites. Thirdly, S-system equations are formed by using the remaining metabolites within the framework of biochemical systems theory. Finally, parameters including rate constants and kinetic orders are estimated by the Levenberg-Marquardt algorithm. The estimation is iterated by setting insignificant kinetic orders at zero, i.e., removing insignificant metabolites. Consequently, a reaction network structure is identified and its mathematical model is obtained. Our approach is validated using a generic inhibition and activation model and its practical application is tested using a simplified model of the glycolysis of Lactococcus lactis MG1363, for which actual time-series data of metabolite concentrations are available. The results indicate the usefulness of our approach and suggest a probable pathway for the production of lactate and acetate. The results also indicate that the approach pinpoints a probable strong inhibition of lactate on the glycolysis pathway.

流通式化学反応器に対する混合拡散モデルを, Lagrange型内挿式を直接微分して導いた式により計算した選点定数を使う修正型の直交選点法によって解き, 反応器の入口, 出口における反応物濃度を基準として数値解の精度を検討した.1次反応における解析解との比較から, Peclet数と反応速度定数が大きな厳しい条件下でも, 単純に選点数を増加させることにより, 倍精度計算において10桁以上の高精度解が得られることがわかった.この結果は, 直交選点法が, Peclet数の値の大きさに応じて適切な刻み幅を選択することが必要な差分法よりも優れていることを明らかに示す.

直交選点法における選点定数を効率よく計算するための方法について検討した.任意の階の選点定数を同一のサブルーチンで計算できるようにするため, それらの再帰性を利用した計算アルゴリズムを提出した.これにより得られた1階, 2階の選点定数の計算値は, 前報の方法による値と同様に高精度を持つことを確認した.また, 三角関数の微分値に対する近似の精度についての検討結果より, 本アルゴリズムが高階の選点定数の計算においてもうまく動作することを確認した.1組の選点定数を計算するのに必要な時間は従来法の時間とほぼ同じであった.

A novel numerical method, based on the combined Taylor series method with a technique for rapidly selecting suitable step sizes, is given for solving two-point boundary value problems. The differential equation for an immobilized enzyme reaction is considered as a model system. Comparisons between numerical and exact solutions for the first-order reaction show that the proposed method promises superhigh-order accuracy that is almost the same as machine accuracy over wide ranges of relevant parameters. The final value of the increment or decrement used to obtain a new estimate in the Newton-Raphson iterative scheme is found to be approximately the same as the relative error of the calculated value, which means that this final value is useful to predict the accuracy of a numerical solution of the differential equation that does not provide an analytical solution. Moreover, it is shown that the selected stepsizes are suitable values that cause little loss-of-significance errors.

Accuracies of the numerical solution of a two-point boundary value problem by the orthogonal collocation method are examined. The collocation constants for first and second derivatives are calculated by the matrix operation method (Method 1) and the method based on Lagrange’s interpolation formula (Method 2). Comparison of the calculated results shows that Method 1 produces marked loss-of-significance errors with increasing number of collocation points, N. These values are used to solve the two-point boundary value problem concerning an immobilized enzyme reaction, and then to calculate the effectiveness factor. These numerical solutions are compared with those of the previously proposed method that provides a numerical solution whose accuracy is almost the same as machine accuracy. It is found that Method 1 provides a numerical solution of significantly low accuracy under the condition of steep concentration gradient and this is never improved, even if N is increased, while Method 2 gives highly accurate values by simply increasing N. These tendencies are shown to depend strongly on the accuracies of the collocation constants.

Numerical differentiation is known to be one of the most difficult numerical calculation methods to obtain reliable calculated values at all times. A simple numerical differentiation method using a combination of finite-difference formulas, derived by approximation of Taylor-series equations, is investigated in order to efficiently perform the sensitivity analysis of large-scale metabolic reaction systems. A result of the application to four basic mathematical functions reveals that the use of the eight-point differentiation formula with a non-dimensionalized stepsize close to 0.01 mostly provides more than 14 digits of accuracy in double precision for the numerical derivatives. Moreover, a result of the application to the modified TCA cycle model indicates that the numerical differentiation method gives the calculated values of steady-state metabolite concentrations within a range of round-off error and also makes it possible to transform the Michaelis-Menten equations into the S-system equations having the kinetic orders whose accuracies are mostly more than 14 significant digits. Because of the simple structure of the numerical differentiation formula and its promising high accuracy, it is evident that the present numerical differentiation method is useful for the analysis of large-scale metabolic reaction systems according to the systematic procedure of BST.

Biochemical systems theory (BST) characterizes a given biochemical system based on the logarithmic gains, rate-constant sensitivities and kinetic-order sensitivities defined at a steady state. This paper describes an efficient method for calculation of the time courses of logarithmic gains, i.e. dynamic logarithmic gains L(Xi, Xj; t), which expresses the percentage change in the value of a dependent variable Xi at a time t in response to an infinitesimal percentage change in the value of an independent variable Xj at t=0. In this method, one first recasts the ordinary differential equations for the dependent variables into an exact canonical nonlinear representation (GMA system) through appropriate transformations of variables. Owing to the structured mathematical form of this representation, the recast system can be fully described by a set of numeric parameters, and the differential equations for the dynamic logarithmic gains can be set up automatically without resource to computer algebra. A simple general-purpose computer program can thus be written that requires only the relevant numeric parameters as input to calculate the time courses of the variables and of the dynamic logarithmic gains for both concentrations and fluxes. Unlike other methods, the proposed method does not require to derive any expression for the partial differentiation of flux expressions with respect to each independent variable. The proposed method has been applied to two kinds of reaction models to elucidate its usefulness.

To determine the usefulness of equations previously proposed for the apparent maximum reaction rate and apparent Michaelis constant of an immobilized enzyme, starch hydrolysis by glucoamylase immobilized on a porous ceramic support was considered as a model system. Initial reaction rates, v0, were measured for a wide range of initial starch concentrations, Sb0, to make a nonlinear plot of Sb0/v0 versus Sb0, and the apparent kinetic parameters were determined from the slopes and intercepts of tangents to the nonlinear plot at given values of Sb0. The equations were found to accurately express the diffusional effect on the kinetic parameters.