Kinetic Ordinary Differential Equations Research Articles

Kinetic Modeling of Cornstalk Cellulose Hydrolysis in Supercritical Water: A Comparative Study of the Effects of Temperature and Residence Time on Derivative Production

Kinetic modeling is essential in understanding and controlling the process of cellulose hydrolysis for producing value-added cellulose derivatives. This study aims to adopt a set of dominate kinetic ordinary differential equations of cornstalk cellulose hydrolysis in supercritical water for mechanism-based prediction of the production of cellulose, glucose, fructose, glyceraldehyde, erythrose, 5-hydroxymethyl furfural, glycolaldehyde, threose, aldose, and other cellulose derivatives from cornstalks under processing conditions with a pressure of 89 MPa and a temperature of 378 °C, as considered in a recent experimental study in the literature. The yield rates of several cellulose derivatives, e.g., glucose, fructose, 5-HMF, and erythrose as predicted by the present model, are close to those of experimental measurements. The model is further used to predict the yield rates of a few new cellulose derivatives, e.g., glycolaldehyde, threose, and aldose, that are potentially generated in cornstalk cellulose hydrolysis in supercritical water. The present model and computational simulations can be utilized as a rational tool to predict, control, and optimize the derivative yields in cellulose hydrolysis in supercritical water via tuning the process parameters, and, therefore, are useful for the optimal production of targeted bio-based fuels and chemicals from cornstalks and other agricultural and municipal wastes.

Accelerating reactive-flow simulations using vectorized chemistry integration

The high cost of chemistry integration is a significant computational bottleneck for realistic reactive-flow simulations using operator splitting. Here we present a methodology to accelerate the solution of the chemical kinetic ordinary differential equations using single-instruction, multiple-data vector processing on CPUs using the OpenCL framework. First, we compared several vectorized integration algorithms using chemical kinetic source terms and analytical Jacobians from the pyJac software against a widely used integration code, CVODEs. Next, we extended the OpenFOAM computational fluid dynamics library to incorporate the vectorized solvers, and we compared the accuracy of a fourth-order linearly implicit integrator–both in vectorized form and a corresponding method native to OpenFOAM—with the community standard chemical kinetics library Cantera. We then applied our methodology to a variety of chemical kinetic models, turbulent intensities, and simulation scales to examine a range of engineering and scientific scale problems, including (pseudo) steady-state as well as time-dependent Reynolds-averaged Navier–Stokes simulations of the Sandia flame D and the Volvo Flygmotor bluff-body stabilized, premixed flame. Subsequently, we compared the performance of the vectorized and native OpenFOAM integrators over the studied models and simulations and found that our vectorized approach performs up to 33–35× faster than the native OpenFOAM solver with high accuracy.

Kinetics-informed neural networks

Chemical kinetics and reaction engineering consists of the phenomenological framework for the disentanglement of reaction mechanisms, optimization of reaction performance and the rational design of chemical processes. Here, we utilize feed-forward artificial neural networks as basis functions to solve ordinary differential equations (ODEs) constrained by differential algebraic equations (DAEs) that describe microkinetic models (MKMs). We present an algebraic framework for the mathematical description and classification of reaction networks, types of elementary reaction, and chemical species. Under this framework, we demonstrate that the simultaneous training of neural nets and kinetic model parameters in a regularized multi-objective optimization setting leads to the solution of the inverse problem through the estimation of kinetic parameters from synthetic experimental data. We analyze a set of scenarios to establish the extent to which kinetic parameters can be retrieved from transient kinetic data, and assess the robustness of the methodology with respect to statistical noise. This approach to inverse kinetic ODEs can assist in the elucidation of reaction mechanisms based on transient data.

Combinatorial mathematical modelling approaches to interrogate rear retraction dynamics in 3D cell migration.

Cell migration in 3D microenvironments is a complex process which depends on the coordinated activity of leading edge protrusive force and rear retraction in a push-pull mechanism. While the potentiation of protrusions has been widely studied, the precise signalling and mechanical events that lead to retraction of the cell rear are much less well understood, particularly in physiological 3D extra-cellular matrix (ECM). We previously discovered that rear retraction in fast moving cells is a highly dynamic process involving the precise spatiotemporal interplay of mechanosensing by caveolae and signalling through RhoA. To further interrogate the dynamics of rear retraction, we have adopted three distinct mathematical modelling approaches here based on (i) Boolean logic, (ii) deterministic kinetic ordinary differential equations (ODEs) and (iii) stochastic simulations. The aims of this multi-faceted approach are twofold: firstly to derive new biological insight into cell rear dynamics via generation of testable hypotheses and predictions; and secondly to compare and contrast the distinct modelling approaches when used to describe the same, relatively under-studied system. Overall, our modelling approaches complement each other, suggesting that such a multi-faceted approach is more informative than methods based on a single modelling technique to interrogate biological systems. Whilst Boolean logic was not able to fully recapitulate the complexity of rear retraction signalling, an ODE model could make plausible population level predictions. Stochastic simulations added a further level of complexity by accurately mimicking previous experimental findings and acting as a single cell simulator. Our approach highlighted the unanticipated role for CDK1 in rear retraction, a prediction we confirmed experimentally. Moreover, our models led to a novel prediction regarding the potential existence of a ‘set point’ in local stiffness gradients that promotes polarisation and rapid rear retraction.

Dynamic optimization of an emulsion copolymerization process for product quality using a deterministic kinetic model with embedded Monte Carlo simulations

We present the dynamic optimization of an emulsion copolymerization process described by a deterministic kinetic ordinary differential equation model including a stochastic Monte Carlo submodel, describing the growth of polymer chains within particles. For the considered semi-batch operation, the time dependent input trajectories for monomer and initiator flow rates are optimized along with the isothermal reactor temperature. We use the surrogate-model-based optimizer MATSuMoTo for the optimization to avoid the need to compute derivatives of the stochastic model. Radial basis functions with linear polynomial tails are selected as surrogate functions which are updated during optimization by the newly evaluated points. A relevant application problem formulation together with results for two case studies are presented. Qualitatively similar input trajectories are obtained for different optimization runs due to the stochastic process and the limited number of iterations. All solutions reduce the batch time significantly.

Synchronization of three electrochemical oscillators: From local to global coupling.

We investigate the formation of synchronization patterns in an oscillatory nickel electrodissolution system in a network obtained by superimposing local and global coupling with three electrodes. We explored the behavior through numerical simulations using kinetic ordinary differential equations, Kuramoto type phase models, and experiments, in which the local to global coupling could be tuned by cross resistances between the three nickel wires. At intermediate coupling strength with predominant global coupling, two of the three oscillators, whose natural frequencies are closer, can synchronize. By adding even a relatively small amount of local coupling (about 9%-25%), a spatially organized partially synchronized state can occur where one of the two synchronized elements is in the center. A formula was derived for predicting the critical coupling strength at which full synchronization will occur independent of the permutation of the natural frequencies of the oscillators over the network. The formula correctly predicts the variation of the critical coupling strength as a function of the global coupling fraction, e.g., with local coupling the critical coupling strength is about twice than that required with global coupling. The results show the importance of the topology of the network on the synchronization properties in a simple three-oscillator setup and could provide guidelines for decrypting coupling topology from identification of synchronization patterns.

Ordinary differential equations and Boolean networks in application to modelling of 6-mercaptopurine metabolism.

We consider two approaches to modelling the cell metabolism of 6-mercaptopurine, one of the important chemotherapy drugs used for treating acute lymphocytic leukaemia: kinetic ordinary differential equations, and Boolean networks supplied with one controlling node, which takes continual values. We analyse their interplay with respect to taking into account ATP concentration as a key parameter of switching between different pathways. It is shown that the Boolean networks, which allow avoiding the complexity of general kinetic modelling, preserve the possibility of reproducing the principal switching mechanism.

Experimental design and inverse problems in plant biological modeling

Abstract. We develop a mathematical and statistical framework to model the actions of underlying metabolites for carbon dioxide assimilation in photosynthesis. To illustrate ideas, we use use a dynamical system model proposed by Zhu, de Sturler and Long (2007) which describes the Calvin Cycle in spinach plants through changes in the concentrations of 38 metabolites (state variables) using non-linear enzyme kinetic ordinary differential equations and mass-balance laws that contain a total of 165 parameters. We address the research question “Once a subset of parameters is fixed and the times at which data is collected are determined, can we identify which metabolites we should measure in order to optimize the confidence in our parameter estimates?” Here, we propose two methods: the ad hoc Simple Linear Regression (SLR) method, and the Optimal Design Criterion method. Both methods suggest promise in determining appropriate sets of metabolites to observe for successful implementation of inverse problems.

A Modified Implicit Iterative Difference Algorithm for Stiff Chemical Kinetic Equations in Complex Combustion System

The splitting-operator method has been widely used in the numerical simulation of complex combustion system with its high computation efficiency. The main difficult encountered in the method is how to integrate the stiff chemical kinetics ordinary differential equation (ODE) efficiently and accurately in each grid node of flow field. Although several ODE software packages such as LSODE and VODE have very predominant performance in integrating stiff and non-stiff case, it couldn’t ensure the positivity of mass fraction in integrating stiff chemical kinetic ODE. A detail liquid rocket hypergolic bipropellant combination of Monomethylhydrazine (MMH)/Red Fuming Nitric Acid (RFNA) chemical mechanism consisting of 550 elementary reactions among 83 species was constituted basing on several relevant literatures firstly. And then a modified implicit iterative difference algorithm with variable time steps was developed to tackle the problem of mass fraction non-positive. The results from the integration of MMH /nitrogen tetroxide (NTO) chemical kinetic ODE indicate that the algorithm is always retain stability and positive value of mass fraction. The time step length would also be automatically adjusted at different chemical reaction time in the algorithm to keep its computational efficiency.

Kinetic modeling and exploratory numerical simulation of chloroplastic starch degradation

BackgroundHigher plants and algae are able to fix atmospheric carbon dioxide through photosynthesis and store this fixed carbon in large quantities as starch, which can be hydrolyzed into sugars serving as feedstock for fermentation to biofuels and precursors. Rational engineering of carbon flow in plant cells requires a greater understanding of how starch breakdown fluxes respond to variations in enzyme concentrations, kinetic parameters, and metabolite concentrations. We have therefore developed and simulated a detailed kinetic ordinary differential equation model of the degradation pathways for starch synthesized in plants and green algae, which to our knowledge is the most complete such model reported to date.ResultsSimulation with 9 internal metabolites and 8 external metabolites, the concentrations of the latter fixed at reasonable biochemical values, leads to a single reference solution showing β-amylase activity to be the rate-limiting step in carbon flow from starch degradation. Additionally, the response coefficients for stromal glucose to the glucose transporter kcat and KM are substantial, whereas those for cytosolic glucose are not, consistent with a kinetic bottleneck due to transport. Response coefficient norms show stromal maltopentaose and cytosolic glucosylated arabinogalactan to be the most and least globally sensitive metabolites, respectively, and β-amylase kcat and KM for starch to be the kinetic parameters with the largest aggregate effect on metabolite concentrations as a whole. The latter kinetic parameters, together with those for glucose transport, have the greatest effect on stromal glucose, which is a precursor for biofuel synthetic pathways. Exploration of the steady-state solution space with respect to concentrations of 6 external metabolites and 8 dynamic metabolite concentrations show that stromal metabolism is strongly coupled to starch levels, and that transport between compartments serves to lower coupling between metabolic subsystems in different compartments.ConclusionsWe find that in the reference steady state, starch cleavage is the most significant determinant of carbon flux, with turnover of oligosaccharides playing a secondary role. Independence of stationary point with respect to initial dynamic variable values confirms a unique stationary point in the phase space of dynamically varying concentrations of the model network. Stromal maltooligosaccharide metabolism was highly coupled to the available starch concentration. From the most highly converged trajectories, distances between unique fixed points of phase spaces show that cytosolic maltose levels depend on the total concentrations of arabinogalactan and glucose present in the cytosol. In addition, cellular compartmentalization serves to dampen much, but not all, of the effects of one subnetwork on another, such that kinetic modeling of single compartments would likely capture most dynamics that are fast on the timescale of the transport reactions.

Analytical and numerical study of photocurrent transients in organic polymer solar cells

This article is an attempt to provide a self consistent picture, including existence analysis and numerical solution algorithms, of the mathematical problems arising from modeling photocurrent transients in organic polymer solar cells (OSCs). The mathematical model for OSCs consists of a system of nonlinear diffusion–reaction partial differential equations (PDEs) with electrostatic convection, coupled to a kinetic ordinary differential equation (ODE). We propose a suitable reformulation of the model that allows us to prove the existence of a solution in both stationary and transient conditions and to better highlight the role of exciton dynamics in determining the device turn-on time. For the numerical treatment of the problem, we carry out a temporal semi-discretization using an implicit adaptive method, and the resulting sequence of differential subproblems is linearized using the Newton–Raphson method with inexact evaluation of the Jacobian. Then, we use exponentially fitted finite elements for the spatial discretization, and we carry out a thorough validation of the computational model by extensively investigating the impact of the model parameters on photocurrent transient times.

Using spline-enhanced ordinary differential equations for PK/PD model development

A spline-enhanced ordinary differential equation (ODE) method is proposed for developing a proper parametric kinetic ODE model and is shown to be a useful approach to PK/PD model development. The new method differs substantially from a previously proposed model development approach using a stochastic differential equation (SDE)-based method. In the SDE-based method, a Gaussian diffusion term is introduced into an ODE to quantify the system noise. In our proposed method, we assume an ODE system with form dx/dt = A(t)x + B(t) where B(t) is a nonparametric function vector that is estimated using penalized splines. B(t) is used to construct a quantitative measure of model uncertainty useful for finding the proper model structure for a given data set. By means of two examples with simulated data, we demonstrate that the spline-enhanced ODE method can provide model diagnostics and serve as a basis for systematic model development similar to the SDE-based method. We compare and highlight the differences between the SDE-based and the spline-enhanced ODE methods of model development. We conclude that the spline-enhanced ODE method can be useful for PK/PD modeling since it is based on a relatively uncomplicated estimation algorithm which can be implemented with readily available software, provides numerically stable, robust estimation for many models, is distribution-free and allows for identification and accommodation of model deficiencies due to model misspecification.

Dynamic reduction of a CH4/air chemical mechanism appropriate for investigating vortex–flame interactions

AbstractThis paper describes two methods, piecewise reusable implementation of solution mapping (PRISM) and dynamic steady‐state approximation (DYSSA), in which chemistry is reduced dynamically to reduce the computational burden in combustion simulations. Each method utilizes the large range in species timescales to reduce the dimensionality to the number of species with slow timescales. The methods are applied within a framework that uses hypercubes to partition multidimensional chemical composition space, where each chemical species concentration, plus temperature, is represented by an axis in space. The dimensionality of the problem is reduced uniquely in each hypercube, but the dimensionality of chemical composition space is not reduced. The dimensionality reduction is dynamic and is different for different hypercubes, thereby escaping the restrictions of global methods in which reductions must be valid for all chemical mixtures. PRISM constructs polynomial equations in each hypercube, replacing the chemical kinetic ordinary differential equation (ODE) system with a set of quadratic polynomials with terms related to the number of species with slow timescales. Earlier versions of PRISM were applied to smaller chemical mechanisms and used all chemical species concentrations as terms. DYSSA is a dynamic treatment of the steady‐state approximation and uses the fast–slow timescale separation to determine the set of steady‐state species in each hypercube. A reduced number of chemical kinetic ODEs are integrated rather than the original full set. PRISM and DYSSA are evaluated for simulations of a pair of counterrotating vortices interacting with a premixed CH4/air laminar flame. DYSSA is sufficiently accurate for use in combustion simulations, and when relative errors are less than 1.0%, speedups on the order of 3 are observed. PRISM does not perform as well as DYSSA with respect to accuracy and efficiency. Although the polynomial evaluation that replaces the ODE solver is sufficiently fast, polynomials are not reused sufficiently to enable their construction cost to be recovered. © 2007 Wiley Periodicals, Inc. 39: 204–220, 2007

Dynamic patterns of gene regulation I: Simple two-gene systems

Regulation of gene activities is studied by means of computer assisted mathematical analysis of ordinary differential equations (ODEs) derived from binding equilibria and chemical reaction kinetics. Here, we present results on cross-regulation of two genes through activator and/or repressor binding. Arbitrary (differentiable) binding function can be used but systematic investigations are presented for gene–regulator complexes with integer valued Hill coefficients up to n = 4 . The dynamics of gene regulation is derived from bifurcation patterns of the underlying systems of kinetic ODEs. In particular, we present analytical expressions for the parameter values at which one-dimensional (transcritical, saddle-node or pitchfork) and/or two-dimensional (Hopf) bifurcations occur. A classification of regulatory states is introduced, which makes use of the sign of a ‘regulatory determinant’ D (being the determinant of the block in the Jacobian matrix that contains the derivatives of the regulator binding functions): (i) systems with D < 0 , observed, for example, if both proteins are activators or repressors, to give rise to one-dimensional bifurcations only and lead to bistability for n ⩾ 2 and (ii) systems with D > 0 , found for combinations of activation and repression, sustain a Hopf bifurcation and undamped oscillations for n > 2 . The influence of basal transcription activity on the bifurcation patterns is described. Binding of multiple subunits can lead to richer dynamics than pure activation or repression states if intermediates between the unbound state and the fully saturated DNA initiate transcription. Then, the regulatory determinant D can adopt both signs, plus and minus.

Computational economy improvements in PRISM

AbstractThe Piecewise Reusable Implementation of Solution Mapping (PRISM) procedure is applied to reactive flow simulations of (9‐species) H2 + air combustion. PRISM takes the solution of the chemical kinetic ordinary differential equation system and parameterizes it with quadratic polynomials. To increase the accuracy, the parameterization is done piecewise, by dividing the multidimensional chemical composition space into hypercubes and constructing polynomials for each hypercube on demand. The polynomial coefficients are stored for subsequent repeated reuse. Initial cost of polynomial construction is expensive, but it recouped as the hypercube is reused, hence computational gain depends on the degree of hypercube reuse. We present two methods that help us to identify hypercubes that will ultimately have high reuse, this being accomplished before the expense of constructing polynomials has been incurred. One method utilizes the rate of movement of the chemical trajectory to estimate the number of steps the trajectory would make through the hypercube. The other method defers polynomial construction until a preset threshold of reuse has been met; an empirical method which, nevertheless, produces a substantial gain. The methods are tested on a 0‐D chemical mixture and reactive flow 1‐D and 2‐D simulations of selected laminar and turbulent H2 + air flames. The computational performance of PRISM is improved by a factor of about 2 for both methods. © 2003 Wiley Periodicals, Inc.*This article is a US Government work and, as such, is in the public domain of the United States of America. Int J Chem Kinet 35: 438–452, 2003

Comparison of different non-linear evaluation methods in thermal analysis

With the aim of determining activation parameters, three computer programs have been investigated. Different routines for the integration of kinetic ordinary differential equations were applied. Sufficient accuracy of the method used is the prerequisite for converging of the algorithm which has to minimize the sum of deviation squares. Our program “TA-kin” enables effective evaluations because of its high flexibility.

TRANSITION WAVES THAT LEAVE BEHIND REGULAR OR IRREGULAR SPATIOTEMPORAL OSCILLATIONS IN A SYSTEM OF THREE REACTION–DIFFUSION EQUATIONS

Transition waves are widespread in the biological and chemical sciences, and have often been successfully modelled using reaction–diffusion systems. I consider a particular system of three reaction–diffusion equations, and I show that transition waves can destabilise as the kinetic ordinary differential equations pass through a Hopf bifurcation, giving rise to either regular or irregular spatiotemporal oscillations behind the advancing transition wave front. In the case of regular oscillations, I show that these are periodic plane waves that are induced by the way in which the transition wave front approaches its terminal steady state. Further, I show that irregular oscillations arise when these periodic plane waves are unstable as reaction–diffusion solutions. The resulting behavior is not related to any chaos in the kinetic ordinary differential equations.