Sensitivity and control analysis of periodically forced reaction networks using the Green's function method

A general sensitivity and control analysis of periodically forced reaction networks with respect to small perturbations in arbitrary network parameters is presented. A well-known property of sensitivity coefficients for periodic processes in dynamical systems is that the coefficients generally become unbounded as time tends to infinity. To circumvent this conceptual obstacle, a relative time or phase variable is introduced so that the periodic sensitivity coefficients can be calculated. By employing the Green's function method, the sensitivity coefficients can be defined using integral control operators that relate small perturbations in the network's parameters and forcing frequency to variations in the metabolite concentrations and reaction fluxes. The properties of such operators do not depend on a particular parameter perturbation and are described by the summation and connectivity relationships within a control-matrix operator equation. The aim of this paper is to derive such a general control-matrix operator equation for periodically forced reaction networks, including metabolic pathways. To illustrate the general method, the two limiting cases of high and low forcing frequency are considered. We also discuss a practically important case where enzyme activities and forcing frequency are modulated simultaneously. We demonstrate the developed framework by calculating the sensitivity and control coefficients for a simple two reaction pathway where enzyme activities enter reaction rates linearly and specifically.     

Crowded surfaces change annealing dynamics of actin filaments

Changes in cell shape that occur in many cellular processes are thought to arise from polymerization of actin filaments near the cell membrane. End-to-end annealing of actin filaments is believed to play only a minor role in this process, as annealing in solution was shown to be a slow process, which is not typical for a bimolecular reaction, its rate constant decreasing over time, being inversely proportional to the filament length. Furthermore, in vitro studies on f-actin solutions were found to display an exponential steady-state length distribution. In the cell, many physiologically important parameters, such as mechanical strength or viscoelastic response are a direct function of the physical properties of the underlying actin cytoskeleton, such as actin filament length distribution and dynamics. How the underlying physical parameters of the actin cytoskeleton may be influenced by the cell surface or molecular crowding remains poorly understood. Using total internal reflection fluorescence (TIRF) microscopy we reinvestigated actin end-to-end annealing in vitro in a more realistic environment. We studied the process near a hydrophilic surface together with crowding agents, in order to mimic the physiological media near the cell membrane, which has substantial amounts of macromolecules present. We find that actin end-to-end annealing changes in three ways near a crowded hydrophilic surface as compared to solution. First the annealing rate becomes a factor of 20 faster than in solution. Second the rate of annealing becomes typical of a bimolecular reaction, shows no length dependence and is basically just a function of the square of the concentration of ends. Lastly the length distribution is Gaussian throughout the entire annealing process. This implicates that dynamic rearrangement of actin filaments by annealing near the leading edge of the cell, could change physical parameters like the mechanical response and contribute significantly to cell motility.     

Confining domains lead to reaction bursts: reaction kinetics in the plasma membrane

Confinement of molecules in specific small volumes and areas within a cell is likely to be a general strategy that is developed during evolution for regulating the interactions and functions of biomolecules. The cellular plasma membrane, which is the outermost membrane that surrounds the entire cell, was considered to be a continuous two-dimensional liquid, but it is becoming clear that it consists of numerous nano-meso-scale domains with various lifetimes, such as raft domains and cytoskeleton-induced compartments, and membrane molecules are dynamically trapped in these domains. In this article, we give a theoretical account on the effects of molecular confinement on reversible bimolecular reactions in a partitioned surface such as the plasma membrane. By performing simulations based on a lattice-based model of diffusion and reaction, we found that in the presence of membrane partitioning, bimolecular reactions that occur in each compartment proceed in bursts during which the reaction rate is sharply and briefly increased even though the asymptotic reaction rate remains the same. We characterized the time between reaction bursts and the burst amplitude as a function of the model parameters, and discussed the biological significance of the reaction bursts in the presence of strong inhibitor activity.     

The evolution of environmental and genetic sex determination in fluctuating environments

Twenty years ago, Bulmer and Bull suggested that disruptive selection, produced by environmental fluctuations, can result in an evolutionary transition from environmental sex determination (ESD) to genetic sex determination (GSD). We investigated the feasibility of such a process, using mutation-limited adaptive dynamics and individual-based computer simulations. Our model describes the evolution of a reaction norm for sex determination in a metapopulation setting with partial migration and variation in an environmental variable both within and between local patches. The reaction norm represents the probability of becoming a female as a function of environmental state and was modeled as a sigmoid function with two parameters, one giving the location (i.e., the value of the environmental variable for which an individual has equal chance of becoming either sex) and the other giving the slope of the reaction norm for that environment. The slope can be interpreted as being set by the level of developmental noise in morph determination, with less noise giving a steeper slope and a more switchlike reaction norm. We found convergence stable reaction norms with intermediate to large amounts of developmental noise for conditions characterized by low migration rates, small differential competitive advantages between the sexes over environments, and little variation between individual environments within patches compared to variation between patches. We also considered reaction norms with the slope parameter constrained to a high value, corresponding to little developmental noise. For these we found evolutionary branching in the location parameter and a transition from ESD toward GSD, analogous to the original analysis by Bulmer and Bull. Further evolutionary change, including dominance evolution, produced a polymorphism acting as a GSD system with heterogamety. Our results point to the role of developmental noise in the evolution of sex determination.     

A modular approach for modeling the cell cycle based on functional response curves

Modeling biochemical reactions by means of differential equations often results in systems with a large number of variables and parameters. As this might complicate the interpretation and generalization of the obtained results, it is often desirable to reduce the complexity of the model. One way to accomplish this is by replacing the detailed reaction mechanisms of certain modules in the model by a mathematical expression that qualitatively describes the dynamical behavior of these modules. Such an approach has been widely adopted for ultrasensitive responses, for which underlying reaction mechanisms are often replaced by a single Hill function. Also time delays are usually accounted for by using an explicit delay in delay differential equations. In contrast, however, S-shaped response curves, which by definition have multiple output values for certain input values and are often encountered in bistable systems, are not easily modeled in such an explicit way. Here, we extend the classical Hill function into a mathematical expression that can be used to describe both ultrasensitive and S-shaped responses. We show how three ubiquitous modules (ultrasensitive responses, S-shaped responses and time delays) can be combined in different configurations and explore the dynamics of these systems. As an example, we apply our strategy to set up a model of the cell cycle consisting of multiple bistable switches, which can incorporate events such as DNA damage and coupling to the circadian clock in a phenomenological way.     

The Relation of Signal Transduction to the Sensitivity and Dynamic Range of Bacterial Chemotaxis

Complex networks of interacting molecular components of living cells are responsible for many important processes, such as signal processing and transduction. An important challenge is to understand how the individual properties of these molecular interactions and biochemical transformations determine the system-level properties of biological functions. Here, we address the issue of the accuracy of signal transduction performed by a bacterial chemotaxis system. The chemotaxis sensitivity of bacteria to a chemoattractant gradient has been measured experimentally from bacterial aggregation in a chemoattractant-containing capillary. The observed precision of the chemotaxis depended on environmental conditions such as the concentration and molecular makeup of the chemoattractant. In a quantitative model, we derived the chemotactic response function, which is essential to describing the signal transduction process involved in bacterial chemotaxis. In the presence of a gradient, an analytical solution is derived that reveals connections between the chemotaxis sensitivity and the characteristics of the signaling system, such as reaction rates. These biochemical parameters are integrated into two system-level parameters: one characterizes the efficiency of gradient sensing, and the other is related to the dynamic range of chemotaxis. Thus, our approach explains how a particular signal transduction property affects the system-level performance of bacterial chemotaxis. We further show that the two parameters can be derived from published experimental data from a capillary assay, which successfully characterizes the performance of bacterial chemotaxis.     

Chiral Polymerization in Open Systems From Chiral-Selective Reaction Rates

We investigate the possibility that prebiotic homochirality can be achieved exclusively through chiral-selective reaction rate parameters without any other explicit mechanism for chiral bias. Specifically, we examine an open network of polymerization reactions, where the reaction rates can have chiral-selective values. The reactions are neither autocatalytic nor do they contain explicit enantiomeric cross-inhibition terms. We are thus investigating how rare a set of chiral-selective reaction rates needs to be in order to generate a reasonable amount of chiral bias. We quantify our results adopting a statistical approach: varying both the mean value and the rms dispersion of the relevant reaction rates, we show that moderate to high levels of chiral excess can be achieved with fairly small chiral bias, below 10%. Considering the various unknowns related to prebiotic chemical networks in early Earth and the dependence of reaction rates to environmental properties such as temperature and pressure variations, we argue that homochirality could have been achieved from moderate amounts of chiral selectivity in the reaction rates.     

Formation of S-nitrosothiols from regiospecific reaction of thiols with N-nitrosotryptophan derivatives

S-nitrosothiols transport nitric oxide in vivo, and so-called transnitrosation reactions (i.e. the transfer of the nitroso function from nitrosothiol to thiolate) are believed to be involved in this process. In the present study we examined the N-nitrosotryptophan derivative-dependent nitrosation of thiols, a hitherto ignored possibility for the formation of S-nitrosothiols. The corresponding products were identified by (15)N-NMR spectrometry. The fact that the reaction proceeded under hypoxic conditions as well as in non-aqueous solution strongly indicated the occurrence of a transnitrosation reaction. Interestingly, S-nitrosothiols could only very slowly transnitrosate N-terminal-blocked tryptophan derivatives like melatonin in non-aqueous solution but did not induce such a reaction in water. The indole moiety of the N-nitrosotryptophan derivatives was fully restituted during the reaction with thiols, as demonstrated by both capillary zone electrophoresis and fluorescence spectroscopy. A determination of the Arrhenius parameters demonstrated that the corresponding rate constants were comparable with the ones known for the transfer of the nitroso function from nitrosothiol to thiolate. Thus, N-nitrosotryptophan-dependent nitrosation of thiols may occur in vivo and might offer the possibility of developing a new class of vasodilative drugs.     

Participation of the Mehler Reaction and Catalase in the Oxygen Exchange of Chloroplast Preparations

We have reinvestigated several aspects of the Mehler reactionin isolated chloroplasts. We have confirmed that the rate ofoxygen uptake is accelerated by a number of biological and artificialelectron carriers. The Mehler reaction also responds to uncouplingagents and exhibits photosynthetic control with ADP. In thepresence of catalase, steady-state oxygen exchange may be establishedand lead to such apparent anomalies as oxygen uptake in thelight followed by oxygen evolution in the dark. The steady stateis a function of catalase concentration, light intensity, andthe presence and concentration of electron carriers, etc., whichaffect the rate of the Mehler reaction.     

Rational polynomial equation as an unbiased approach for the kinetic studies of Drosophila melanogaster acetylcholinesterase reaction mechanism

The hydrolysis of substrates by cholinesterases does not follow the Michaelis–Menten reaction mechanism. The well-known inhibition by excess substrate is often accompanied by an unexpectedly high activity at low substrate concentrations. It appears that these peculiarities are the consequence of an unusual architecture of the active site, which conducts the substrate molecule over many stages before it is cleaved and released. Structural and kinetic data also suggest that two substrate molecules can attach at the same time to the free, as well as to the acetylated, enzyme. We present a procedure which provides an unbiased framework for mathematical modelling of such complex reaction mechanisms. It is based on regression analysis of a rational polynomial using classical initial rate data. The determination of polynomial degree reveals the number of independent parameters that can be evaluated from the available information. Once determined, these parameters can substantially facilitate the construction and evaluation of a kinetic model reflecting the expected molecular events in an enzymic reaction. We also present practical suggestions for testing the postulated kinetic model, using an original thermodynamic approach and an isolated effect in a specifically mutated enzyme.     

Multimerization-cyclization of DNA fragments as a method of conformational analysis

Ligation of short DNA fragments results in the formation of linear and circular multimers of various lengths. The distribution of products in such a reaction is often used to evaluate fragment bending caused by specific chemical modification, by bound ligands or by the presence of irregular structural elements. We have developed a more rigorous quantitative approach to the analysis of such experimental data based on determination of j-factors for different multimers from the distribution of the reaction products. j-Factors define the effective concentration of one end of a linear chain in the vicinity of the other end. To extract j-factors we assumed that kinetics of the reaction is described by a system of differential equations where j-factors appear as coefficients. The assumption was confirmed by comparison with experimental data obtained here for DNA fragments containing A-tracts. At the second step of the analysis j-factors are used to determine conformational parameters of DNA fragments: the equilibrium bend angle, the bending rigidity of the fragment axis, and the total twist of the fragments. This procedure is based on empirical equations that connect the conformational parameters with the set of j-factors. To obtain the equations, we computed j-factors for a large array of conformational parameters that describe model fragments. The approach was tested on both simulated and actual experimental data for DNA fragments containing A-tracts. A-tract DNA bend angle determined here is in good agreement with previously published data. We have established a set of experimental conditions necessary for the data analysis to be successful.     

Aldehyde PEGylation kinetics: a standard protein versus a pharmaceutically relevant single chain variable fragment

The mPEG-aldehyde PEGylation with two different PEG sizes and two proteins was experimentally determined with respect to yield, conversion, and selectivity. The kinetic behavior of these PEGylation reactions was simulated using a numerically solved set of differential equations. We show that the assumption of an inactivation of mPEG-aldehyde is crucial for the simulation of the overall PEGylation and that the inactivation is pH-dependent. We further demonstrate that ideal PEGylation parameters such as pH, temperature, reaction time, and protein concentration need to be chosen carefully depending on the protein and PEG size. In terms of selectivity and yield, we show that the reaction should be stopped before the highest mono-PEG concentration is reached. Moreover, room temperature and a slightly acidic pH of approximately 6 are good starting points. In conclusion, selectivity can be optimized choosing a shorter reaction time and a reduced reaction temperature.     

Trunk orientation causes asymmetries in leg function in small bird terrestrial locomotion

In contrast to the upright trunk in humans, trunk orientation in most birds is almost horizontal (pronograde). It is conceivable that the orientation of the heavy trunk strongly influences the dynamics of bipedal terrestrial locomotion. Here, we analyse for the first time the effects of a pronograde trunk orientation on leg function and stability during bipedal locomotion. For this, we first inferred the leg function and trunk control strategy applied by a generalized small bird during terrestrial locomotion by analysing synchronously recorded kinematic (three-dimensional X-ray videography) and kinetic (three-dimensional force measurement) quail locomotion data. Then, by simulating quail gaits using a simplistic bioinspired numerical model which made use of parameters obtained in in vivo experiments with real quail, we show that the observed asymmetric leg function (left-skewed ground reaction force and longer leg at touchdown than at lift-off) is necessary for pronograde steady-state locomotion. In addition, steady-state locomotion becomes stable for specific morphological parameters. For quail-like parameters, the most common stable solution is grounded running, a gait preferred by quail and most of the other small birds. We hypothesize that stability of bipedal locomotion is a functional demand that, depending on trunk orientation and centre of mass location, constrains basic hind limb morphology and function, such as leg length, leg stiffness and leg damping.     

Optimal experiment selection for parameter estimation in biological differential equation models

ABSTRACT: BACKGROUND: Parameter estimation in biological models is a common yet challenging problem. In this work we explore the problem for gene regulatory networks modeled by differential equations with unknown parameters, such as decay rates, reaction rates, Michaelis-Menten constants, and Hill coefficients. We explore the question to what extent parameters can be efficiently estimated by appropriate experimental selection. RESULTS: A minimization formulation is used to find the parameter values that best fit the experiment data. When the data is insufficient, the minimization problem often has many local minima that fit the data reasonably well. We show that selecting a new experiment based on the local Fisher Information of one local minimum generates additional data that allows one to successfully discriminate among the many local minima. The parameters can be estimated to high accuracy by iteratively performing minimization and experiment selection. We show that the experiment choices are roughly independent of which local minima is used to calculate the local Fisher Information. CONCLUSIONS: We show that by an appropriate choice of experiments, one can, in principle, efficiently and accurately estimate all the parameters of gene regulatory network. In addition, we demonstrate that appropriate experiment selection can also allow one to restrict model predictions without constraining the parameters using many fewer experiments. We suggest that predicting model behaviors and inferring parameters represent two different approaches to model calibration with different requirements on data and experimental cost.     

Comparison of Two New Robust Parameter Estimation Methods for the Power Function Distribution

Estimation of any probability distribution parameters is vital because imprecise and biased estimates can be misleading. In this study, we investigate a flexible power function distribution and introduced new two methods such as, probability weighted moments, and generalized probability weighted methods for its parameters. We compare their results with L-moments, trimmed L-moments by a simulation study and a real data example based on performance measures such as, mean square error and total deviation. We concluded that all the methods perform well in the case of large sample size (n>30), however, the generalized probability weighted moment method performs better for small sample size.     

A novel cost function to estimate parameters of oscillatory biochemical systems

Oscillatory pathways are among the most important classes of biochemical systems with examples ranging from circadian rhythms and cell cycle maintenance. Mathematical modeling of these highly interconnected biochemical networks is needed to meet numerous objectives such as investigating, predicting and controlling the dynamics of these systems. Identifying the kinetic rate parameters is essential for fully modeling these and other biological processes. These kinetic parameters, however, are not usually available from measurements and most of them have to be estimated by parameter fitting techniques. One of the issues with estimating kinetic parameters in oscillatory systems is the irregularities in the least square (LS) cost function surface used to estimate these parameters, which is caused by the periodicity of the measurements. These irregularities result in numerous local minima, which limit the performance of even some of the most robust global optimization algorithms. We proposed a parameter estimation framework to address these issues that integrates temporal information with periodic information embedded in the measurements used to estimate these parameters. This periodic information is used to build a proposed cost function with better surface properties leading to fewer local minima and better performance of global optimization algorithms. We verified for three oscillatory biochemical systems that our proposed cost function results in an increased ability to estimate accurate kinetic parameters as compared to the traditional LS cost function. We combine this cost function with an improved noise removal approach that leverages periodic characteristics embedded in the measurements to effectively reduce noise. The results provide strong evidence on the efficacy of this noise removal approach over the previous commonly used wavelet hard-thresholding noise removal methods. This proposed optimization framework results in more accurate kinetic parameters that will eventually lead to biochemical models that are more precise, predictable, and controllable.     

Poly(silicate)-metabolizing silicatein in siliceous spicules and silicasomes of demosponges comprises dual enzymatic activities (silica polymerase and silica esterase)

Siliceous sponges can synthesize poly(silicate) for their spicules enzymatically using silicatein. We found that silicatein exists in silica-filled cell organelles (silicasomes) that transport the enzyme to the spicules. We show for the first time that recombinant silicatein acts as a silica polymerase and also as a silica esterase. The enzymatic polymerization/polycondensation of silicic acid follows a distinct course. In addition, we show that silicatein cleaves the ester-like bond in bis(p-aminophenoxy)-dimethylsilane. Enzymatic parameters for silica esterase activity are given. The reaction is completely blocked by sodium hexafluorosilicate and E-64. We consider that the dual function of silicatein (silica polymerase and silica esterase) will be useful for the rational synthesis of structured new silica biomaterials.     

The contribution of proximity and orientation to catalytic reaction rates

We present calculations of the possible magnitude of propinquity which has been proposed to play an important role in enzymic catalysis. The effect has been evaluated as in the past by calculating the ratio of bimolecular to intramolecular reaction rates. The ratio is estimated for intramolecular catalysis in rigid systems as well as for systems with five and six rotatable bonds. Our method differs from others mainly in the way cyclization has been treated. The reaction rate (in all cases) is proportional to the probability that the reactive units have the appropriate spatial and orientational positioning for reaction. This probability is obtained by evaluating a distance distribution function within the spatial and angular intervals to which the units are constrained after reaction. For the bimolecular case we have made the usual assumption that the distribution function is uniform. For the intramolecular reaction, neither the spatial nor the angular part of the distribution function is uniform. The pertinent parameters in this case are the bond lengths and angles, and the statistical weight matrices describing torsional rotation. The difficulty in obtaining analytical expressions for the distribution function is circumvented by using Monte Carlo methods. It is argued that the spatial contribution to rate accelerations in rigid systems can be as high as 10⁷ M, depending upon the size of the volume to which the reactive units are constrained after reaction. The limitation on the smallest physically reasonable volume is estimated from considerations of energy requirements and vibrational amplitudes. Accelerations by five- and six-membered ring cyclizations were estimated at 10³ M, the six-membered ring exhibiting the smaller rate enhancement.     

A graphical simulation model of the entire DNA process associated with the analysis of short tandem repeat loci

The use of expert systems to interpret short tandem repeat DNA profiles in forensic, medical and ancient DNA applications is becoming increasingly prevalent as high-throughput analytical systems generate large amounts of data that are time-consuming to process. With special reference to low copy number (LCN) applications, we use a graphical model to simulate stochastic variation associated with the entire DNA process starting with extraction of sample, followed by the processing associated with the preparation of a PCR reaction mixture and PCR itself. Each part of the process is modelled with input efficiency parameters. Then, the key output parameters that define the characteristics of a DNA profile are derived, namely heterozygote balance (Hb) and the probability of allelic drop-out p(D). The model can be used to estimate the unknown efficiency parameters, such as π_extraction. ‘What-if’ scenarios can be used to improve and optimize the entire process, e.g. by increasing the aliquot forwarded to PCR, the improvement expected to a given DNA profile can be reliably predicted. We demonstrate that Hb and drop-out are mainly a function of stochastic effect of pre-PCR molecular selection. Whole genome amplification is unlikely to give any benefit over conventional PCR for LCN.