Parameter optimization by using differential elimination: a general approach for introducing constraints into objective functions

Background

The investigation of network dynamics is a major issue in systems and synthetic biology. One of the essential steps in a dynamics investigation is the parameter estimation in the model that expresses biological phenomena. Indeed, various techniques for parameter optimization have been devised and implemented in both free and commercial software. While the computational time for parameter estimation has been greatly reduced, due to improvements in calculation algorithms and the advent of high performance computers, the accuracy of parameter estimation has not been addressed.

Results

We propose a new approach for parameter optimization by using differential elimination, to estimate kinetic parameter values with a high degree of accuracy. First, we utilize differential elimination, which is an algebraic approach for rewriting a system of differential equations into another equivalent system, to derive the constraints between kinetic parameters from differential equations. Second, we estimate the kinetic parameters introducing these constraints into an objective function, in addition to the error function of the square difference between the measured and estimated data, in the standard parameter optimization method. To evaluate the ability of our method, we performed a simulation study by using the objective function with and without the newly developed constraints: the parameters in two models of linear and non-linear equations, under the assumption that only one molecule in each model can be measured, were estimated by using a genetic algorithm (GA) and particle swarm optimization (PSO). As a result, the introduction of new constraints was dramatically effective: the GA and PSO with new constraints could successfully estimate the kinetic parameters in the simulated models, with a high degree of accuracy, while the conventional GA and PSO methods without them frequently failed.

Conclusions

The introduction of new constraints in an objective function by using differential elimination resulted in the drastic improvement of the estimation accuracy in parameter optimization methods. The performance of our approach was illustrated by simulations of the parameter optimization for two models of linear and non-linear equations, which included unmeasured molecules, by two types of optimization techniques. As a result, our method is a promising development in parameter optimization.

     

A fast model of voltage-dependent NMDA receptors

NMDA receptors are among the crucial elements of central nervous system models. Recent studies show that both conductance and kinetics of these receptors are changing voltage-dependently in some parts of the brain. Therefore, several models have been introduced to simulate their current. However, on the one hand, kinetic models—which are able to simulate these voltage-dependent phenomena—are computationally expensive for modeling of large neural networks. On the other hand, classic exponential models, which are computationally less expensive, are not able to simulate the voltage-dependency of these receptors, accurately. In this study, we have modified these classic models to endow them with the voltage-dependent conductance and time constants. Temperature sensitivity and desensitization of these receptors are also taken into account. We show that, it is possible to simulate the most important physiological aspects of NMDA receptor’s behavior using only three to four differential equations, which is significantly smaller than the previous kinetic models. Consequently, it seems that our model is both fast and physiologically plausible and therefore is a suitable candidate for the modeling of large neural networks.     

Kinetic Model Development for Accelerated Stability Studies

Accelerated stability coupled with modeling to predict the stability of compounds, blends, and products at long-term storage conditions provides significant benefits in science-based decision-making throughout drug substance and drug product development. The study can often be completed, including data analysis in the space of three working weeks, and the information gathered and learning made in this time period can rival years of traditional analysis. The speed of the studies allows an earlier assessment of risk to quality enabling appropriate risk mitigation strategies to be implemented in a timely manner. The scientific foundation is based upon Arrhenius kinetic equations that can be linear or nonlinear in time, and can be based upon water vapor pressure or liquid water activity (relative humidity). A variety of kinetic models are evaluated, and the best model is chosen based upon both Bayesian information criteria and an automated assessment of kinetic model parameters fitting within acceptable ranges. Confidence intervals are estimated based upon a bootstrapping approach. Moisture vapor transmission rate models are applied on top of the resulting kinetic models in order to simulate different packaging types and the use of desiccant. The kinetic models are integrated with the prediction of packaging humidity over time to create a long-term prediction of impurities and other phenomena. The resulting models have been shown to be useful for not only the prediction of drug product impurities in long-term storage but other physical phenomena as well such as hydrate development and solvate loss.     

Kinetic model of Proteus mirabilis swarm colony development

 Proteus mirabilis colonies display striking symmetry and periodicity. Based on experimental observations of cellular differentiation and group motility, a kinetic model has been developed to describe the swarmer cell differentiation-dedifferentiation cycle and the spatial evolution of swimmer and swarmer cells during Proteus mirabilis swarm colony development. A key element of the model is the age dependence of swarmer cell behaviour, in particular specifying a minimal age for motility and maximum age for septation and dedifferentiation to swimmer cells. Density thresholds for collective motility by mature swarmer cells serve to synchronize the movements of distinct swarmer cell groups and thus help provide temporal coherence to colony expansion cycles. Numerical computations show that the model fits experimental data by generating a complete swarming plus consolidation cycle period that is robust to changes in parameters which affect other aspects of swarmer cell migration and colony development. The kinetic equations underlying this model provide a different mathematical basis for a temporal oscillator from reaction-diffusion partial differential equations. The modelling shows that Proteus colony geometries arise as a consequence of macroscopic rules governing collective motility. Thus, in this case, pattern formation results from the operation of an adaptive bacterial system for spreading on solid substrates, not as an independent biological function. Kinetic models similar to this one may be applicable to periodic phenomena displayed by other biological systems with differentiated components of defined lifetimes. Received 3 July 1996; received in revised form 9 December 1996     

Multiscale bio-chemo-mechanical model of intimal hyperplasia

We consider a computational multiscale framework of a bio-chemo-mechanical model for intimal hyperplasia. With respect to existing models, we investigate the interactions between hemodynamics, cellular dynamics and biochemistry on the development of the pathology. Within the arterial wall, we propose a mathematical model consisting of kinetic differential equations for key vascular cell types, collagen and growth factors. The luminal hemodynamics is modeled with the Navier–Stokes equations. Coupling hypothesis among time and space scales are proposed to build a tractable modeling of such a complex multifactorial and multiscale pathology. A one-dimensional numerical test-case is presented for validation by comparing the results of the framework with experiments at short and long timescales. Our model permits to capture many cellular phenomena which have a central role in the physiopathology of intimal hyperplasia. Results are quantitatively and qualitatively consistent with experimental findings at both short and long timescales.

     

Modelling atypical CYP3A4 kinetics: principles and pragmatism

The Michaelis-Menten model, and the existence of a single active site for the interaction of substrate with drug metabolizing enzyme, adequately describes a substantial number of in vitro metabolite kinetic data sets for both clearance and inhibition determination. However, in an increasing number of cases (involving most notably, but not exclusively, CYP3A4), atypical kinetic features are observed, e.g., auto- and heteroactivation; partial, cooperative, and substrate inhibition; concentration-dependent effector responses (activation/inhibition); limited substrate substitution and inhibitory reciprocity necessitating sub-group classification. The phenomena listed above cannot be readily interpreted using single active site models and the literature indicates that three types of approaches have been adopted. First the 'nai ve' approach of using the Michaelis-Menten model regardless of the kinetic behaviour, second the 'empirical' approach (e.g., employing the Hill or uncompetitive inhibition equations to model homotropic phenomena of sigmoidicity and substrate inhibition, respectively) and finally, the 'mechanistic' approach. The later includes multisite kinetic models derived using the same rapid equilibrium/steady-state assumptions as the single-site model. These models indicate that 2 or 3 binding sites exist for a given CYP3A4 substrate and/or effector. Multisite kinetic models share common features, depending on the substrate kinetics and the nature of the effector response observed in vitro, which allow a generic model to be proposed. Thus although more complex than the other two approaches, they show more utility and can be comprehensively applied in relatively simple versions that can be readily generated from generic model. Multisite kinetic features, observed in isolated hepatocytes as well as in microsomes from hepatic tissue and heterologous expression systems, may be evident in substrate depletion-time profiles as well as in metabolite formation rates. Failure to adequately account for multisite kinetic phenomena will compromise any attempts to predict human drug clearance and drug-drug interaction potential from in vitro data.     

Phenomenological and molecular-level Petri net modeling and simulation of long-term potentiation

Petri net-based modeling methods have been used in many research projects to represent biological systems. Among these, the hybrid functional Petri net (HFPN) was developed especially for biological modeling in order to provide biologists with a more intuitive Petri net-based method. In the literature, HFPNs are used to represent kinetic models at the molecular level. We present two models of long-term potentiation previously represented by differential equations which we have transformed into HFPN models: a phenomenological synapse model and a molecular-level model of the CaMKII regulation pathway. Through simulation, we obtained results similar to those of previous studies using these models. Our results open the way to a new type of modeling for systems biology where HFPNs are used to combine different levels of abstraction within one model. This approach can be useful in fully modeling a system at the molecular level when kinetic data is missing or when a full study of a system at the molecular level it is not within the scope of the research.     

Modeling and artificial intelligence approaches to enzyme systems

Modeling is a means of formulating and testing complex hypotheses. Useful modeling is now possible with biological laboratory microcomputers with which experimenters feel comfortable. Artificial intelligence (AI) is sufficiently similar to modeling that AI techniques, now becoming usable on microcomputers, are applicable to modeling. Microcomputer and AI applications to physiological system studies with multienzyme models and with kinetic models of isolated enzymes are described. Using an IBM PC microcomputer, we have been able to fit kinetic enzyme models; to extend this process to design kinetic experiments by determining the optimal conditions; and to construct an enzyme (hexokinase) kinetics data base. We have also used a PC to do most of the constructing of complex multienzyme models, initially with small simple BASIC programs; alternative methods with standard spreadsheet or data base programs have been defined. Formulating and solving differential equations in appropriate representational languages, and sensitivity analysis, are soon likely to be feasible with PCs. Much of the modeling process can be stated in terms of AI expert systems, using sets of rules for fitting and evaluating models and designing further experiments. AI techniques also permit critiquing and evaluating the data, experiments, and hypotheses being modeled, and can be extended to supervise the calculations involved.     

Models of spatially restricted biochemical reaction systems

Many reactions within the cell occur only in specific intracellular regions. Such local reaction networks give rise to microdomains of activated signaling components. The dynamics of microdomains can be visualized by live cell imaging. Computational models using partial differential equations provide mechanistic insights into the interacting factors that control microdomain dynamics. The mathematical models show that, for membrane-initiated signaling, the ratio of the surface area of the plasma membrane to the volume of the cytoplasm, the topology of the signaling network, the negative regulators, and kinetic properties of key components together define microdomain dynamics. Thus, patterns of locally restricted signaling reaction systems can be considered an emergent property of the cell.     

A comparative investigation of certain difference equations and related differential equations: Implications for model-building

Many mathematical models for physical and biological problems have been and will be built in the form of differential equations or systems of such equations. With the advent of digital computers one has been able to find (approximate) solutions for equations that used to be intractable. Many of the mathematical techniques used in this area amount to replacing the given differential equations by appropriate difference equations, so that extensive research has been done into how to choose appropriate difference equations whose solutions are “good” approximations to the solutions of the given differential equations. The present paper investigates a different, although related problem. For many physical and biological phenomena the “continuum” type of thinking, that is at the basis of any differential equation, is not natural to the phenomenon, but rather constitutes an approximation to a basically discrete situation: in much work of this type the “infinitesimal step lengths” handled in the reasoning which leads up to the differential equation, are not really thought of as infinitesimally small, but as finite; yet, in the last stage of such reasoning, where the differential equation rises from the differentials, these “infinitesimal” step lengths are allowed to go to zero: that is where the above-mentioned approximation comes in. Under this kind of circumstances, it seems more natural tobuild themodel as adiscrete difference equation (recurrence relation) from the start, without going through the painful, doubly approximative process of first, during the modeling stage, finding a differential equation to approximate a basically discrete situation, and then, for numerical computing purposes, approximating that differential equation by a difference scheme. The paper pursues this idea for some simple examples, where the old differential equation, though approximative in principle, had been at least qualitatively successful in describing certain phenomena, and shows that this idea, though plausible and sound in itself, does encounter some difficulties. The reason is that each differential equation, as it is set up in the way familiar to theoretical physicists and biologists, does correspond to a plethora of discrete difference equations, all of which in the limit (as step length→0) yield the same differential equation, but whose solutions, for not too small step length, are often widely different, some of them being quite irregular. The disturbing thing is that all these difference equations seem to adequately represent the same (physical or biological) reasoning as the differential equation in question. So, in order to choose the “right” difference equation, one may need to draw upon more detailed (physical or) biological considerations. All this does not say that one should not prefer discrete models for phenomena that seem to call for them; but only that their pursuit may require additional (physical or) biological refinement and insight. The paper also investigates some mathematical problems related to the fact of many difference equations being associated with one differential equation.     

Practical identifiability of growth and substrate consumption models

The estimation of parameters in several dynamic models, which describe growth and substrate consumption, has been carried out using a modified Gauss-Newton-type method. The four models considered are Monod, Contois, linear specific growth rate, and an enzyme kinetic model. The initial values of the differential equations are included in the parameter vector which will be estimated. The efficiency of the method and the confidence limits of the parameters were studied using simulated measurement noise. The experimental results describe Trichoderma viride growing on glucose as the main carbon source.     

Kinetic model for dynamic response of three-phase fluidized bed biofilm reactor for wastewater treatment

Step changes in inlet concentration has been introduced into the completely mixed three-phase fluidized bed biofilm reactor treating simulated domestic wastewater to study the dynamic behavior of the system and to establish the suitable kinetic model from the response curve. Three identical reactors having different biomass volumes were operated in parallel. It was found that the response curves showed second-order characteristics, and thus at least two first-order differential equations are necessary to simulate the substrate and biomass response curves. Nonlinear regression analysis was performed using different types of rate equations and their corresponding kinetic parameters were used to simulate the theoretical response curve using the Runge–Kutta numerical integration method. As a result, although various types of conventional biokinetic models such as Monod, Haldane and Andrew types were examined, all the theoretical substrate response curves underestimated time constants compared to the actual substrate response plots. On the other hand, the theoretical curve of the kinetic model that incorporates adsorption term has best fit to the actual response in most of the cases. Thus, it was concluded that adsorption of substrate onto biofilm and carrier particles has significant effect on the dynamic response in biofilm processes.     

A critical review of mathematical models and data used in diabetology

The literature dealing with mathematical modelling for diabetes is abundant. During the last decades, a variety of models have been devoted to different aspects of diabetes, including glucose and insulin dynamics, management and complications prevention, cost and cost-effectiveness of strategies and epidemiology of diabetes in general. Several reviews are published regularly on mathematical models used for specific aspects of diabetes. In the present paper we propose a global overview of mathematical models dealing with many aspects of diabetes and using various tools. The review includes, side by side, models which are simple and/or comprehensive; deterministic and/or stochastic; continuous and/or discrete; using ordinary differential equations, partial differential equations, optimal control theory, integral equations, matrix analysis and computer algorithms.     

From molecular networks to qualitative cell behavior

Adaptation and behavior are characteristics of life which are fundamentally dynamic. If we want to model the living cell we have to describe it as a dynamic system. Typical dynamic models are based on quantitative differential equations requiring very detailed kinetic knowledge. Alternative modeling techniques for less fine-grained information are better suited to available functional genomics data. As such, constraint-based techniques and qualitative modeling have proven themselves to be valid approaches in cell biology. These approaches offer formal support to check the consistency of molecular networks against phenotypic observations in the light of dynamic systems.     

Mathematical models of threshold phenomena in the nerve membrane

The types of mathematical model which have been used to represent all-or-none behavior in the nerve membrane may be classified as follows: (1) thediscontinuous threshold phenomenon, in which differential equations with discontinuous functions provide both a discontinuity of response as a function of stimulus intensity at threshold and a finite maximum latency, (2) thesingular-point threshold phenomenon which exists in a phase space having analytic functions in its differential equations and having a singular point with one characteristic root positive and the rest with negative real parts, the latency being unbounded, and (3) thequasi threshold phenomenon, which has a finite maximum latency and continuous functions, but neither a true discontinuity in response nor an exact threshold. Several models of the nerve membrane in the literature are classified accordingly, and the applicability of the different types of threshold phenomena to the membrane is discussed, including an extension to a stochastic model.     

A Bond Graph Model of the Cardiovascular System

The study of the autonomic nervous system (ANS) function has shown to provide useful indicators for risk stratification and early detection on a variety of cardiovascular pathologies. However, data gathered during different tests of the ANS are difficult to analyse, mainly due to the complex mechanisms involved in the autonomic regulation of the cardiovascular system (CVS). Although model-based analysis of ANS data has been already proposed as a way to cope with this complexity, only a few models coupling the main elements involved have been presented in the literature. In this paper, a new model of the CVS, representing the ventricles, the circulatory system and the regulation of the CVS activity by the ANS, is presented. The models of the vascular system and the ventricular activity have been developed using the Bond Graph formalism, as it proposes a unified representation for all energetic domains, facilitating the integration of mechanic and hydraulic phenomena. In order to take into account the electro-mechanical behaviour of both ventricles, an electrophysiologic model of the cardiac action potential, represented by a set of ordinary differential equations, has been integrated. The short-term ANS regulation of heart rate, cardiac contractility and peripheral vasoconstriction is represented by means of continuous transfer functions. These models, represented in different continuous formalisms, are coupled by using a multi-formalism simulation library. Results are presented for two different autonomic tests, namely the Tilt Test and the Valsalva Manoeuvre, by comparing real and simulated signals.     

Mathematical models of hierarchically structured cell populations under equilibrium with application to the epidermis

Abstract. There are three categories of keratinocytes in the germinative compartment of the epidermis – stem, transit-amplifying and post-mitotic. Their population structure is hierarchical. This means that stem cells differentiate into transit-amplifying cells which, after a few rounds of division, become post-mitotic cells. The cell processes of birth, differentiation, death and migration affect the composition and proliferation rate of the germinative compartment. These phenomena are quantified by various cell kinetic parameters. In this paper we derive equations that relate these parameters for different models of hierarchically structured cell populations in equilibrium. We include in the models asymmetric and symmetric division, variations in cell-cycle times, apoptosis and variation in the number of transit generations. We conclude that variation in cell-cycle times need only be considered if apoptosis is not negligible. If it is negligible, then only average cell-cycle times are needed. Unfortunately, it is impossible to predict the importance of apoptosis from the available experimental data. However, the strength of its effect is determined by the other parameters, especially the fraction of cycling stem cells. We show that variation in the number of transit generations can have a potentially large effect on cell birth rate. We also show that cell birth rate does not directly depend on the mean transit-amplifying cell-cycle time, only on the mean stem cell-cycle time. We argue that 'homogeneous cell population' equations should not be used to study hierarchical cell populations as has been done in the past. Finally we argue that stem cell parameters and transit-amplifying cell parameters should not be lumped together.     

Kinetic studies of enzymatic hydrolysis of insoluble cellulose: Derivation of a mechanistic kinetic model

A comprehensive mechanistic kinetic model for enzymatic hydrolysis of insoluble cellulose has been synthesized by combining models for several key aspects which have been derived independent of each other. The model takes into account the major contributing factors: the nature of the enzyme system, the structure of cellulose, and the mode of interaction between the enzyme and cellulose molecules. It consists of a set of simultaneously occurring ordinary differential equations with ten kinetic constants. All of the kinetic constants have been determined independently by carrying out critically designed experiments, and they appear in the comprehensive model without any arbitrary manipulations. The governing equations of the model have been numerically simulated by means of the computer subroutine CSMP III. The model predicts the progress of hydrolysis of cellulose over a wide range of experimental conditions and hydrolysis times reasonably well. The model can even be applied to predict the progress of hydrolysis for intensively pretreated cellulose with a minor adjustment. The applicability of the model for the actual process development is also discussed.