Highly accurate computation of dynamic sensitivities in metabolic reaction systems by a Taylor series method

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.     

A new method of identifying a systems based on the minimal quadratic discrepancy criterion for biophysics problems

A wide range of biophysical systems are described by nonlinear dynamic models mathematically presented as a set of ordinary differential equations in the Cauchy explicit form: [formula: see text] Fij(X1(t),..,XN(t),t), (i = 1,...,N, j = 1,...,M), where Fij (X1(t), ..., XN(t), t) is a set of basis functions satisfying the Lipschitz condition. We investigate the problem of evaluation of model constants aij (the system identification) using experimental data about the time dependence of the dynamic parameters of the system Xi(t). A new method of system identification for the class of similar nonlinear dynamic models is proposed. It is shown that the problem of identifying an initial nonlinear model can be reduced to the solution of a system of linear equations for the matrix of the dynamic model constants [aj]i. It is proposed to determine the set of dynamic model constants aij using the criterion of minimal quadratic discrepancy for the time dependence of the set of dynamic parameters Xi(t). An important special case of the nonlinear model, the quadratic model, is considered. Test problems of identification using this method are presented for two nonlinear systems: the Van der Pol type multiparametric nonlinear oscillator and the strange attractor of Ressler, a widely known example of dynamic systems showing the stochastic behavior.     

崔-Lawson和Logistic方程参数的优化估计方法

本文提出用单形寻优与微分方程数值解法的联合方法,进行生态学中一些微分动力系统的参数的优化估计。用这种方法来估计崔-Lawson和Logistic方程的各参数效果极好。     

崔—Lawson和Logistic方程参数的优化估计方法

本文提出用单形寻优与微分方程数值解法的联合方法,进行生态学中一些微分动力系统的参数的优化估计。用这种方法来估计崔-Lawson和Logistic方程的各参数效果极好。     

Three algorithms for interpreting models consisting of ordinary differential equations: Sensitivity coefficients,sensitivity functions,global optimization

This paper reports the development and application of three powerful algorithms for the analysis and simulation of mathematical models consisting of ordinary differential equations. First, we describe an extended parameter sensitivity analysis: we measure the relative sensitivities of many dynamical behaviors of the model to perturbations of each parameter. We check sensitivities to parameter variation over both small and large ranges. These two extensions of a common technique have applications in parameter estimation and in experimental design. Second, we compute sensitivity functions, using an efficient algorithm requiring just one model simulation to obtain all sensitivities of state variables to all parameters as functions of time. We extend the analysis to a behavior which is not a state variable. Third, we present an unconstrained global optimization algorithm, and apply it in a novel way: we determine the input to the model, given an optimality criterion and typical outputs. The algorithm itself is an efficient one for high-order problems, and does not get stuck at local extrema. We apply the sensitivity analysis, sensitivity functions, and optimization algorithm to a sixth-order nonlinear ordinary differential equation model for human eye movements. This application shows that the algorithms are not only practicable for high-order models, but also useful as conceptual tools.     

A quantitative model of myosin phosphorylation and the photomechanical response of the isolated sphincter pupillae of the frog iris.

The time courses of isometrically recorded photomechanical responses of isolated sphincter pupillae of Rana pipiens can be accurately predicted by a set of differential equations derived from phosphorylation theory of smooth muscle contraction. We compared actual light-stimulated contractions with calculated ones over a wide range of stimulus intensities (56-fold) and durations (0.4-4.0 s). The hypothetical Ca++-calmodulin-myosin light chain kinase cascade acts as a "valve" to control the flow of ATP through a phosphorylation-dephosphorylation cycle. When the rate of flow of ATP through the phosphorylation-dephosphorylation cycle is increased, the percentage of phosphorylated myosin increases. The time courses of the concentrations of phosphorylated myosin during different responses are seen to be functions of the time courses of the opening and closing of the coupling cascade "valve." The calculations predict experimentally measurable intermediate variables, which can aid the investigation of the application of quantitative phosphorylation theory to amphibian sphincter pupillae and to smooth muscle in general.     

Perception of effort in isometric and dynamic muscular contraction.

The perception of muscular effort was studied using estimation and production methods in the adductor pollicis and quadriceps. A psychometric scale (percentage magnitude) was used. Static contractions were studied in the adductor pollicis, and both dynamic (isokinetic) and static contractions were studied in quadriceps. Linear and logarithmic equations were fitted for the perceived effort as a percentage of the maximum in relation to the produced percentage maximal force or torque. The logarithmic exponent was around or above 1.0. No significant difference was found between mean exponent and intercept values for the adductor pollicis and the quadriceps, or when estimated or produced values for the two muscles were compared. There was no difference in the same subjects between the equations for static and dynamic contractions with low angular velocity of the quadriceps.     

Parameter estimation in Hodgkin-Huxley-Type equations for membrane action potentials in nerve and heart muscle

A method is presented which renders parameter estimation possible in systems of non-linear differential equations where normally no solution exists in terms of analytic functions and which have to be solved numerically. The method uses the concept of sensitivity equations. Two examples are given, taking mathematical models for membrane action potentials in nerve and heart muscle by Hodgkin and Huxley and by Beeler and Reuter. The model equations together with the corresponding system of sensitivity equations are given, which are necessary to estimate maximum conductivity coefficients defining the interactions of different ionic current components. A computer program is described and results of action potential numerical analysis are presented using simulated data. It can be seen, that even with superimposed simulated noise the real parameter values are estimated in an excellent manner. The method can be used to interpret observed changes in action potential time courses under physiological and pharmacological conditions.     

PLMaddon: a power-law module for the Matlab SBToolbox

PLMaddon is a General Public License (GPL) software module designed to expand the current version of the SBToolbox (a Matlab toolbox for systems biology; www.sbtoolbox.org) with a set of functions for the analysis of power-law models, a specific class of kinetic models, set in ordinary differential equations (ODE) and in which the kinetic orders can have positive/negative non-integer values. The module includes functions to generate power-law Taylor expansions of other ODE models (e.g. Michaelis-Menten type models), as well as algorithms to estimate steady-states. The robustness and sensitivity of the models can also be analysed and visualized by computing the power-law's logarithmic gains and sensitivities.     

Rational engineering of the TOL meta-cleavage pathway

The meta-cleavage pathway of Pseudomonas putida mt-2 was simulated using a biochemical systems simulation developed by Regan (1996). A non-competitive inhibition term for catechol-2,3-dioxygenase (C23O) by 2-OH-pent-2,4-dienoate (Ki = 150 μM) was incorporated into the model. The simulation predicted steady state accumulation levels in the μM range for metabolites pre-meta-cleavage, and in the mM range for metabolites post-meta-cleavage. The logarithmic gains L[V-i, Xj] and L[X-i, Xj] clearly indicated that the pathway was most sensitive to the concentration of the starting substrate, benzoate, and the first enzyme of the pathway, toluate-1, 2-dioxygenase (TO). The simulation was validated experimentally; it was found that the amplification of TO increased the steady state flux from 0.024 to 0.091 (mmol/g cell dwt)/h. This resulted in an increased accumulation of a number of the pathway metabolites (intra- and extracellularly), especially cis-diol, 4-OH-2-oxovalerate, and 4-oxalocrotonate. Metabolic control analysis indicated that C23O was, in fact, the major controling enzymic step of the pathway with a scaled control coefficient of 0.83. The amplification of TO resulted in a shift of some of the control away from C23O. Catechol-2,3-dioxygenase, however, remained as the major controling element of the pathway. Copyright 1998 John Wiley & Sons, Inc.     

A reliable Taylor series-based computational method for the calculation of dynamic sensitivities in large-scale metabolic reaction systems: Algorithm and software evaluation

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.     

Calculation errors of time-varying flux control coefficients obtained from elasticity coefficients by means of summation and connectivity theorems in metabolic control analysis

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.     

A variational approach to parameter estimation in ordinary differential equations

ABSTRACT: BACKGROUND: Ordinary differential equations are widely-used in the field of systems biology andchemical engineering to model chemical reaction networks. Numerous techniques havebeen developed to estimate parameters like rate constants, initial conditions or steady stateconcentrations from time-resolved data. In contrast to this countable set of parameters, theestimation of entire courses of network components corresponds to an innumerable set ofparameters. RESULTS: The approach presented in this work is able to deal with course estimation for extrinsicsystem inputs or intrinsic reactants, both not being constrained by the reaction networkitself. Our method is based on variational calculus which is carried out analytically toderive an augmented system of differential equations including the unconstrainedcomponents as ordinary state variables. Finally, conventional parameter estimation isapplied to the augmented system resulting in a combined estimation of courses andparameters. CONCLUSIONS: The combined estimation approach takes the uncertainty in input courses correctly intoaccount. This leads to precise parameter estimates and correct confidence intervals. Inparticular this implies that small motifs of large reaction networks can be analysedindependently of the rest. By the use of variational methods, elements from control theoryand statistics are combined allowing for future transfer of methods between the two fields.     

Canonical sensitivities: a useful tool to deal with large perturbations in metabolic network modeling

The dynamic range of metabolic models can be extended to deal with large perturbations by introducing the related concepts of "generalized" kinetic order and "canonical" sensitivities. Generalized kinetic orders are built as a well-defined non linear combination of the canonical sensitivities coefficients, which in turn are obtained by a least-squares regression on central composite factorial design data. In a such way, the whole domain of the operating variables is mapped without need to determine locally neither the first nor the second order model derivatives. The method was validated through numerical simulations, its predictions being compared with those coming from a Michaelis-Menten formalism taken as reference. In parallel, two variants of the Power-law formalism (S-system, least-squares GMA) also were tested. The canonical sensitivities method produced the widest range to predict metabolite concentrations and metabolic fluxes at the steady states. In addition, the variation pattern for the logarithmic gains and for the characteristic eigenvalues have been accurately determined from a unique overall model, being both required to make realistic analysis in metabolic engineering. The achieved information also can be expressed in terms of those typical coefficients derived from the Metabolic Control Analysis (MCA). Even if current first order Power-law or MCA formalisms were used, the canonical sensitivities approach provides a significant advantage, since complete sets of homologous, accurate, locally valid metabolic coefficients can be simultaneously recovered from the array proposed, being representative of the whole range of the operating variables instead of a unique nominal condition as is usual.     

Simplified dynamics model of planar two-joint arm movements

A theoretical framework is presented that describes a way in which the inverse dynamics equations of motion of planar two-joint arm movements (EX-model) are reformulated in a simple form. A single point was assumed to define both the wrist and elbow joint centers, and thus the motion of two points in extrinsic space was represented by second-order differential equations to provide the variables in the reformulation (RE-) model. Through an analytical processes, it was shown that the RE-model for reproducing the shoulder joint torque consists of the linearly scaled moment per unit mass responsible for accelerating the wrist and elbow points about the shoulder joint, while that for reproducing the elbow joint torque consists of the linearly scaled moment per unit mass responsible for accelerating the wrist point about the elbow. The scaling factors for variables in the RE-model were based solely on the values for segment lengths, while in the EX-model the inertial parameter data for the segments are involved in its representation. The inertial parameter data of six-arm specimens from the cadaver experiment of Chandler et al. (1975, AMRL Technical Report, Wright-Patterson Air Force Base, OH) were used to develop and verify the numeric solutions of the RE-model. The adequacy of the model varied somewhat among subjects, but minor changes of the physical parameters of the arm segments enabled perfect reformulation, regardless of the specimens. The potential abilities of the RE-model to deal with the complexities in motor control with more simple control schemes are discussed.     

Artificial neural network as an alternative to multiple regression analysis in optimizing formulation parameters of cytarabine liposomes

The objective of the study was to optimize the formulation parameters of cytarabine liposomes by using artificial neural networks (ANN) and multiple regression analysis using 3(3) factorial design (FD). As model formulations, 27 formulations were prepared. The formulation variables, drug (cytarabine)/lipid (phosphatidyl choline [PC] and cholesterol [Chol]) molar ratio (X1), PC/Chol in percentage ratio of total lipids (X2), and the volume of hydration medium (X3) were selected as the independent variables; and the percentage drug entrapment (PDE) was selected as the dependent variable. A set of causal factors was used as tutorial data for ANN and fed into a computer. The optimization was performed by minimizing the generalized distance between the predicted values of each response and the optimized one that was obtained individually. In case of 3(3) factorial design, a second-order full-model polynomial equation and a reduced model were established by subjecting the transformed values of independent variables to multiple regression analysis, and contour plots were drawn using the equation. The optimization methods developed by both ANN and FD were validated by preparing another 5 liposomal formulations. The predetermined PDE and the experimental data were compared with predicted data by paired t test, no statistically significant difference was observed. ANN showed less error compared with multiple regression analysis. These findings demonstrate that ANN provides more accurate prediction and is quite useful in the optimization of pharmaceutical formulations when compared with the multiple regression analysis method.     

Smooth bistable S-systems

S-systems have been used as models of biochemical systems for over 30 years. One of their hallmarks is that, although they are highly non-linear, their steady states are characterised by linear equations. This allows streamlined analyses of stability, sensitivities and gains as well as objective, mathematically controlled comparisons of similar model designs. Regular S-systems have a unique steady state at which none of the system variables is zero. This makes it difficult to represent switching phenomena, as they occur, for instance, in the expression of genes, cell cycle phenomena and signal transduction. Previously, two strategies were proposed to account for switches. One was based on a technique called recasting, which permits the modelling of any differentiable non-linearities, including bistability, but typically does not allow steady-state analyses based on linear equations. The second strategy formulated the switching system in a piece-wise fashion, where each piece consisted of a regular S-system. A representation gleaned from a simplified form of recasting is proposed and it is possible to divide the characterisation of the steady states into two phases, the first of which is linear, whereas the other is non-linear, but easy to execute. The article discusses a representative pathway with two stable states and one unstable state. The pathway model exhibits strong separation between the stable states as well as hysteresis.     

The Use of a Symbol Processing Computer Language in Point Estimation of Parameters and in Construction of Confidence Intervals

This paper presents a method to generate automatically computer programs which are necessary for parameter estimation, hypothesis tests and construction of confidence intervals by the maximum likelihood method. The spectral or density function of the random variable is arbitrary, but must be known and given in closed form. The programming language used is the symbol processing language LIBAFORM, whose statements are interpreted by a package of LISP-routines. The application of the method is illustrated by the analysis of a linear model whose residuals follow a logarithmic F-distribution, and the analysis of a dose-response curve.     

An efficient simulation method for discrete-value controlled large-scale neuromyoskeletal system models

An efficient Euler-Adams hybrid integration scheme for simulating on the computer discrete-value controlled large-scale neuromyoskeletal system models is presented. If, as discussed in the model, the differential equations describing the recruitment and excitation dynamics of the muscular subsystem are independent of the corresponding contraction-dynamical state variables, they can be integrated separately over certain time intervals by a modified Euler routine that handles discontinuous right-hand sides efficiently. The resulting myostates can then be stored and used as continuous input values for the subsequent integration by an Adams predictor-corrector algorithm of the remaining contraction-dynamical and skeletomechanical state differential equations. With such an Euler-Adams hybrid integration routine one avoids the detrimental effects and efficiency losses associated with frequent stop-restart cycles of otherwise efficient Adams-type algorithms, which cycles are forced by discontinuities on the right-hand side of the myostate equations. In the example presented, a reduction in the execution time by a factor of about 5 could be achieved by implementing the proposed technique.