Numerical bifurcation analysis of delay differential equations arising from physiological modeling

This paper has a dual purpose. First, we describe numerical methods for continuation and bifurcation analysis of steady state solutions and periodic solutions of systems of delay differential equations with an arbitrary number of fixed, discrete delays. Second, we demonstrate how these methods can be used to obtain insight into complex biological regulatory systems in which interactions occur with time delays: for this, we consider a system of two equations for the plasma glucose and insulin concentrations in a diabetic patient subject to a system of external assistance. The model has two delays: the technological delay of the external system, and the physiological delay of the patient's liver. We compute stability of the steady state solution as a function of two parameters, compare with analytical results and compute several branches of periodic solutions and their stability. These numerical results allow to infer two categories of diabetic patients for which the external system has different efficiency.     

Numerical bifurcation analysis of the pattern formation in a cell based auxin transport model

Transport models of growth hormones can be used to reproduce the hormone accumulations that occur in plant organs. Mostly, these accumulation patterns are calculated using time step methods, even though only the resulting steady state patterns of the model are of interest. We examine the steady state solutions of the hormone transport model of Smith et al. (Proc Natl Acad Sci USA 103(5):1301–1306, 2006) for a one-dimensional row of plant cells. We search for the steady state solutions as a function of three of the model parameters by using numerical continuation methods and bifurcation analysis. These methods are more adequate for solving steady state problems than time step methods. We discuss a trivial solution where the concentrations of hormones are equal in all cells and examine its stability region. We identify two generic bifurcation scenarios through which the trivial solution loses its stability. The trivial solution becomes either a steady state pattern with regular spaced peaks or a pattern where the concentration is periodic in time.     

Temperature and burn injury prediction of human skin exposed to microwaves: a model analysis

A one-dimensional multi-layer model is presented to characterize skin temperature rises and burn processes resulting from skin exposure to microwaves. Temperature variations and damage function analyses in the skin tissue exposed to microwaves were predicted depending on blood perfusion rate, thermal conductivity, power density, and exposure time. Thermal wave model was applied and the bio-heat transfer equation was solved using the finite difference time domain method. The thermal wave model of bio-heat transfer predicts a lower temperature rise than a model that uses Pennes’ equation. When approaching steady state, the solutions overlaps with that obtained using the Pennes’ equation. The results obtained may help to analyze the consequences of short-time high-power MW exposures in biological tissues.     

Transfer of Natural Micro Structures to Bionic Lightweight Design Proposals

The abstraction of complex biological lightweight structure features into a producible technical component is a funda- mental step within the transfer of design principles from nature to technical lightweight solutions. A major obstacle for the transfer of natural lightweight structures to technical solutions is their peculiar geometry. Since natural lightweight structures possess irregularities and often have extremely complex forms due to elaborate growth processes, it is usually necessary to simplify their design principles. This step of simplification/abstraction has been used in different biomimetic methods, but so far, it has an arbitrary component, i.e. it crucially depends on the competence of the person who executes the abstraction. This paper describes a new method for abstraction and specialization of natural micro structures for technical lightweight compo- nents. The new method generates stable lightweight design principles by using topology optimization within a design space of preselected biological archetypes such as diatoms or radiolarian. The resulting solutions are adapted to the technical load cases and production processes, can be created in a large variety, and may be further optimized e.g. by using parametric optimization.     

Multistationarity in Sequential Distributed Multisite Phosphorylation Networks

Multisite phosphorylation networks are encountered in many intracellular processes like signal transduction, cell-cycle control, or nuclear signal integration. In this contribution, networks describing the phosphorylation and dephosphorylation of a protein at n sites in a sequential distributive mechanism are considered. Multistationarity (i.e., the existence of at least two positive steady state solutions of the associated polynomial dynamical system) has been analyzed and established in several contributions. It is, for example, known that there exist values for the rate constants where multistationarity occurs. However, nothing else is known about these rate constants. Here, we present a sign condition that is necessary and sufficient for multistationarity in n-site sequential, distributive phosphorylation. We express this sign condition in terms of linear systems, and show that solutions of these systems define rate constants where multistationarity is possible. We then present, for n≥2, a collection of feasible linear systems, and hence give a new and independent proof that multistationarity is possible for n≥2. Moreover, our results allow to explicitly obtain values for the rate constants where multistationarity is possible. Hence, we believe that, for the first time, a systematic exploration of the region in parameter space where multistationarity occurs has become possible. One consequence of our work is that, for any pair of steady states, the ratio of the steady state concentrations of kinase-substrate complexes equals that of phosphatase-substrate complexes.     

Axes of Extrapolation in Risk Assessment

Extrapolation in risk assessment involves the use of data and information to estimate or predict something that has not been measured or observed. Reasons for extrapolation include that the number of combinations of environmental stressors and possible receptors is too large to characterize risks comprehensively, that direct characterization is sometimes impossible, and that the power to characterize risk in a particular situation can be enhanced by using information obtained in other similar situations. Three types of extrapolation are common in risk assessments: biological (including between taxa and across levels of biological organization), temporal, and spatial. They can be thought of conceptually as the axes of a 3-dimensional graph defining the state space of biological, temporal, and spatial scales within which extrapolations are made. Each of these types of extrapolation can introduce uncertainties into risk assessments. Such uncertainties may be reduced through synergistic research facilitated by the sharing of methods, models, and data used by human health and ecological scientists     

The influence of link protein stabilization on the viscometric properties of proteoglycan aggregate solutions

The dynamic, steady-shear and transient shear flow properties of precisely prepared link-stable (s0 136, 66% aggregate) and link-free (s0 93, 59% aggregate) proteoglycan aggregate solutions at concentrations ranging from 10 to 50 mg/ml were determined using a cone-on-plate viscometer in a mechanical spectrometer. All proteoglycan solutions tested possessed: (1) linear viscoelastic properties - as measured by the dynamic complex modulus under small amplitude steady oscillatory conditions (1 less than or equal to omega less than or equal to 100 rad/s) - and (2) nonlinear shear-rate dependent apparent viscosities and primary normal stress difference under steady shearing conditions (0.25 less than or equal to gamma less than or equal to 250 s-1). Our transient flow data show that all proteoglycan aggregate solutions exhibited transient stress overshoot effects in shear stress and normal stress. From these steady and transient flow data, we conclude that link protein stabilized aggregates have significant effects on their dynamic and steady-shear properties as well as transient flow properties. The transient stress overshoot data provide a measure of the energy per unit volume of fluid required to overcome the proteoglycan networks in solution from a resting state. Thus we found that link-stable aggregates form much stronger networks than link-free aggregates. This is corroborated by the fact that link-stable aggregates form more elastic (lower than delta) and stiffer (higher [G*]) networks than link-free aggregates. The complete spectrum of viscometric flow data is entirely compatible with the proposed role of link protein in adding structural stability to the proteoglycan-hyaluronate bond. In cartilage, the enhanced strength of the networks formed by link-stable aggregates may play an important role in determining the material properties of the tissue and thereby contribute to the functional capacity of cartilage in diarthrodial joints.     

Repeated small perturbation approach reveals transcriptomic steady states

The study of biological systems dynamics requires elucidation of the transitions of steady states. A "small perturbation" approach can provide important information on the "steady state" of a biological system. In our experiments, small perturbations were generated by applying a series of repeating small doses of ultraviolet radiation to a human keratinocyte cell line, HaCaT. The biological response was assessed by monitoring the gene expression profiles using cDNA microarrays. Repeated small doses (10 J/m2) of ultraviolet B (UVB) exposure modulated the expression profiles of two groups of genes in opposite directions. The genes that were up-regulated have functions mainly associated with anti-proliferation/anti-mitogenesis/apoptosis, and the genes that were down-regulated were mainly related to proliferation/mitogenesis/anti-apoptosis. For both groups of genes, repetition of the small doses of UVB caused an immediate response followed by relaxation between successive small perturbations. This cyclic pattern was suppressed when large doses (233 or 582.5 J/m2) of UVB were applied. Our method and results contribute to a foundation for computational systems biology, which implicitly uses the concept of steady state.     

Stability of transformation systems representable by tree graphs: Extension to structures of biological importance

It is known that systems representable by tree graphs have entirely real eigenvalues near a steady state (Hyver, 1973). Here it is shown that the eigenvalues are negative, thus ensuring local stability.The method used in the proof allows some extensions which may be of considerable biological importance in certain cases, for example where linear systems containing circuits are involved, or enzymatic reactions, or autocatalytic reactions.     

An Efficient, Nonlinear Stability Analysis for Detecting Pattern Formation in Reaction Diffusion Systems

Reaction diffusion systems are often used to study pattern formation in biological systems. However, most methods for understanding their behavior are challenging and can rarely be applied to complex systems common in biological applications. I present a relatively simple and efficient, nonlinear stability technique that greatly aids such analysis when rates of diffusion are substantially different. This technique reduces a system of reaction diffusion equations to a system of ordinary differential equations tracking the evolution of a large amplitude, spatially localized perturbation of a homogeneous steady state. Stability properties of this system, determined using standard bifurcation techniques and software, describe both linear and nonlinear patterning regimes of the reaction diffusion system. I describe the class of systems this method can be applied to and demonstrate its application. Analysis of Schnakenberg and substrate inhibition models is performed to demonstrate the methods capabilities in simplified settings and show that even these simple models have nonlinear patterning regimes not previously detected. The real power of this technique, however, is its simplicity and applicability to larger complex systems where other nonlinear methods become intractable. This is demonstrated through analysis of a chemotaxis regulatory network comprised of interacting proteins and phospholipids. In each case, predictions of this method are verified against results of numerical simulation, linear stability, asymptotic, and/or full PDE bifurcation analyses.     

Controllability of non-linear biochemical systems

Mathematical methods of biochemical pathway analysis are rapidly maturing to a point where it is possible to provide objective rationale for the natural design of metabolic systems and where it is becoming feasible to manipulate these systems based on model predictions, for instance, with the goal of optimizing the yield of a desired microbial product. So far, theory-based metabolic optimization techniques have mostly been applied to steady-state conditions or the minimization of transition time, using either linear stoichiometric models or fully kinetic models within biochemical systems theory (BST). This article addresses the related problem of controllability, where the task is to steer a non-linear biochemical system, within a given time period, from an initial state to some target state, which may or may not be a steady state. For this purpose, BST models in S-system form are transformed into affine non-linear control systems, which are subjected to an exact feedback linearization that permits controllability through independent variables. The method is exemplified with a small glycolytic-glycogenolytic pathway that had been analyzed previously by several other authors in different contexts.     

Recurrent model of the dependence of a recessive lethal mutation on fitness components and its solution

Experimental modeling and subsequent analysis of microevolutionary processes often involves estimation of fitness components, such as male mating competitiveness, female fecundity, progeny viability, meiotic drive upon the formation of sexual products, and so on. In this connection, we have developed a new method for estimating the relationship among fitness components and their effect on selection on frequencies of the mutant individuals in several generations, which can be applied to the case of rapid elimination of a lethal mutation from the population. The method of estimating unknown fitness components is based on the known estimates of other components and a relationship among the known and unknown components calculated from the frequency dynamics. Using the new method, a biological explanation of non-uniqueness of the admissible solutions. The method employs analysis of the form of the admissible solution region (at preset confidence intervals for the observed frequencies) in the space of the possible values.     

Peculiarities of state diagrams of aqueous solutions of cryoprotective agents

The phase transitions in aqueous solutions of glycerol and PEO-1500 within the temperature range of +30 to −150 °C have been studied using the methods of thermoplastic analysis and volumetric scanning tensodilatometry. We present the revealed phenomenon of cluster cyrystallization of these solutions as well as principles of describing this phenomenon using state diagrams, containing the intervals of concentration corresponding to the existence of amorphous and cryocolloid fractions. We note that for the cryocolloid fraction, a low temperature association of molecules of cryoprotective agents leads the formation of ice nanocrystals either close to or directly inside the aggregations. These fractions exist in cooled cryoprotective solutions near the vitrification temperatures of the liquid phase and may contribute to the initiation of damaging events in cryopreserved biological systems. Our data may be helpful in explaining the peculiarities observed during crystallization of cryoprotective solutions and may further contribute to a broader understanding of the principles of protection and protocol optimization of biological materials at temperatures approaching vitrification.     

The steady state kinetics of some biological systems: I

The steady state kinetics of some typical catalytic systems of biological importance have been formulated. The conditions for the existence of a maximum or limiting velocity are examined and discussed. In particular it is shown that the limiting velocity for a given component is simply the rate expression for a given number of steps of theoverall process; from the general condition for a limiting velocity these steps may be specified. The stringency of the conditions which must be imposed upon the steady state solution in order that it may be assumed that one or more steps are essentially at equilibrium is pointed out. The application of the general method to coupled or branched systems and to cyclic systems is briefly discussed. This work was done while the author was a Senior Research Fellow of the National Institute of Health.     

Efficient Characterization of Parametric Uncertainty of Complex (Bio)chemical Networks

Parametric uncertainty is a particularly challenging and relevant aspect of systems analysis in domains such as systems biology where, both for inference and for assessing prediction uncertainties, it is essential to characterize the system behavior globally in the parameter space. However, current methods based on local approximations or on Monte-Carlo sampling cope only insufficiently with high-dimensional parameter spaces associated with complex network models. Here, we propose an alternative deterministic methodology that relies on sparse polynomial approximations. We propose a deterministic computational interpolation scheme which identifies most significant expansion coefficients adaptively. We present its performance in kinetic model equations from computational systems biology with several hundred parameters and state variables, leading to numerical approximations of the parametric solution on the entire parameter space. The scheme is based on adaptive Smolyak interpolation of the parametric solution at judiciously and adaptively chosen points in parameter space. As Monte-Carlo sampling, it is “non-intrusive” and well-suited for massively parallel implementation, but affords higher convergence rates. This opens up new avenues for large-scale dynamic network analysis by enabling scaling for many applications, including parameter estimation, uncertainty quantification, and systems design.     

A computational framework for the design of optimal protein synthesis

Despite the establishment of design principles to optimize codon choice for heterologous expression vector design, the relationship between codon sequence and final protein yield remains poorly understood. In this work, we present a computational framework for the identification of a set of mutant codon sequences for optimized heterologous protein production, which uses a codon-sequence mechanistic model of protein synthesis. Through a sensitivity analysis on the optimal steady state configuration of protein synthesis we are able to identify the set of codons, that are the most rate limiting with respect to steady state protein synthesis rate, and we replace them with synonymous codons recognized by charged tRNAs more efficient for translation, so that the resulting codon-elongation rate is higher. Repeating this procedure, we iteratively optimize the codon sequence for higher protein synthesis rate taking into account multiple constraints of various types. We determine a small set of optimized synonymous codon sequences that are very close to each other in sequence space, but they have an impact on properties such as ribosomal utilization or secondary structure. This limited number of sequences can then be offered for further experimental study. Overall, the proposed method is very valuable in understanding the effects of the different properties of mRNA sequences on the final protein yield in heterologous protein production and it can find applications in synthetic biology and biotechnology.     

Systems biology towards life in silico: mathematics of the control of living cells

Systems Biology is the science that aims to understand how biological function absent from macromolecules in isolation, arises when they are components of their system. Dedicated to the memory of Reinhart Heinrich, this paper discusses the origin and evolution of the new part of systems biology that relates to metabolic and signal-transduction pathways and extends mathematical biology so as to address postgenomic experimental reality. Various approaches to modeling the dynamics generated by metabolic and signal-transduction pathways are compared. The silicon cell approach aims to describe the intracellular network of interest precisely, by numerically integrating the precise rate equations that characterize the ways macromolecules’ interact with each other. The non-equilibrium thermodynamic or ‘lin–log’ approach approximates the enzyme rate equations in terms of linear functions of the logarithms of the concentrations. Biochemical Systems Analysis approximates in terms of power laws. Importantly all these approaches link system behavior to molecular interaction properties. The latter two do this less precisely but enable analytical solutions. By limiting the questions asked, to optimal flux patterns, or to control of fluxes and concentrations around the (patho)physiological state, Flux Balance Analysis and Metabolic/Hierarchical Control Analysis again enable analytical solutions. Both the silicon cell approach and Metabolic/Hierarchical Control Analysis are able to highlight where and how system function derives from molecular interactions. The latter approach has also discovered a set of fundamental principles underlying the control of biological systems. The new law that relates concentration control to control by time is illustrated for an important signal transduction pathway, i.e. nuclear hormone receptor signaling such as relevant to bone formation. It is envisaged that there is much more Mathematical Biology to be discovered in the area between molecules and Life.     

Quantifying dynamic stability of genetic memory circuits

Bistability/Multistability has been found in many biological systems including genetic memory circuits. Proper characterization of system stability helps to understand biological functions and has potential applications in fields such as synthetic biology. Existing methods of analyzing bistability are either qualitative or in a static way. Assuming the circuit is in a steady state, the latter can only reveal the susceptibility of the stability to injected DC noises. However, this can be inappropriate and inadequate as dynamics are crucial for many biological networks. In this paper, we quantitatively characterize the dynamic stability of a genetic conditional memory circuit by developing new dynamic noise margin (DNM) concepts and associated algorithms based on system theory. Taking into account the duration of the noisy perturbation, the DNMs are more general cases of their static counterparts. Using our techniques, we analyze the noise immunity of the memory circuit and derive insights on dynamic hold and write operations. Considering cell-to-cell variations, our parametric analysis reveals that the dynamic stability of the memory circuit has significantly varying sensitivities to underlying biochemical reactions attributable to differences in structure, time scales, and nonlinear interactions between reactions. With proper extensions, our techniques are broadly applicable to other multistable biological systems.