Mathematical identification of critical reactions in the interlocked feedback model

Dynamic simulations are necessary for understanding the mechanism of how biochemical networks generate robust properties to environmental stresses or genetic changes. Sensitivity analysis allows the linking of robustness to network structure. However, it yields only local properties regarding a particular choice of plausible parameter values, because it is hard to know the exact parameter values in vivo. Global and firm results are needed that do not depend on particular parameter values. We propose mathematical analysis for robustness (MAR) that consists of the novel evolutionary search that explores all possible solution vectors of kinetic parameters satisfying the target dynamics and robustness analysis. New criteria, parameter spectrum width and the variability of solution vectors for parameters, are introduced to determine whether the search is exhaustive. In robustness analysis, in addition to single parameter sensitivity analysis, robustness to multiple parameter perturbation is defined. Combining the sensitivity analysis and the robustness analysis to multiple parameter perturbation enables identifying critical reactions. Use of MAR clearly identified the critical reactions responsible for determining the circadian cycle in the Drosophila interlocked circadian clock model. In highly robust models, while the parameter vectors are greatly varied, the critical reactions with a high sensitivity are uniquely determined. Interestingly, not only the per-tim loop but also the dclk-cyc loop strongly affect the period of PER, although the dclk-cyc loop hardly changes its amplitude and it is not potentially influential. In conclusion, MAR is a powerful method to explore wide parameter space without human-biases and to link a robust property to network architectures without knowing the exact parameter values. MAR identifies the reactions critically responsible for determining the period and amplitude in the interlocked feedback model and suggests that the circadian clock intensively evolves or designs the kinetic parameters so that it creates a highly robust cycle.     

Exploring the gap between dynamic and constraint-based models of metabolism

Systems biology provides new approaches for metabolic engineering through the development of models and methods for simulation and optimization of microbial metabolism. Here we explore the relationship between two modeling frameworks in common use namely, dynamic models with kinetic rate laws and constraint-based flux models. We compare and analyze dynamic and constraint-based formulations of the same model of the central carbon metabolism of Escherichia coli. Our results show that, if unconstrained, the space of steady states described by both formulations is the same. However, the imposition of parameter-range constraints can be mapped into kinetically feasible regions of the solution space for the dynamic formulation that is not readily transferable to the constraint-based formulation. Therefore, with partial kinetic parameter knowledge, dynamic models can be used to generate constraints that reduce the solution space below that identified by constraint-based models, eliminating infeasible solutions and increasing the accuracy of simulation and optimization methods.     

A data-informed mean-field approach to mapping of cortical parameter landscapes

Constraining the many biological parameters that govern cortical dynamics is computationally and conceptually difficult because of the curse of dimensionality. This paper addresses these challenges by proposing (1) a novel data-informed mean-field (MF) approach to efficiently map the parameter space of network models; and (2) an organizing principle for studying parameter space that enables the extraction biologically meaningful relations from this high-dimensional data. We illustrate these ideas using a large-scale network model of the Macaque primary visual cortex. Of the 10-20 model parameters, we identify 7 that are especially poorly constrained, and use the MF algorithm in (1) to discover the firing rate contours in this 7D parameter cube. Defining a “biologically plausible” region to consist of parameters that exhibit spontaneous Excitatory and Inhibitory firing rates compatible with experimental values, we find that this region is a slightly thickened codimension-1 submanifold. An implication of this finding is that while plausible regimes depend sensitively on parameters, they are also robust and flexible provided one compensates appropriately when parameters are varied. Our organizing principle for conceptualizing parameter dependence is to focus on certain 2D parameter planes that govern lateral inhibition: Intersecting these planes with the biologically plausible region leads to very simple geometric structures which, when suitably scaled, have a universal character independent of where the intersections are taken. In addition to elucidating the geometry of the plausible region, this invariance suggests useful approximate scaling relations. Our study offers, for the first time, a complete characterization of the set of all biologically plausible parameters for a detailed cortical model, which has been out of reach due to the high dimensionality of parameter space.     

Feedback-mediated dynamics in a model of a compliant thick ascending limb

The tubuloglomerular feedback (TGF) system in the kidney, which is a key regulator of filtration rate, has been shown in physiologic experiments in rats to mediate oscillations in tubular fluid pressure and flow, and in NaCl concentration in the tubular fluid of the thick ascending limb (TAL). In this study, we developed a mathematical model of the TGF system that represents NaCl transport along a TAL with compliant walls. The model was used to investigate the dynamic behaviors of the TGF system. A bifurcation analysis of the TGF model equations was performed by deriving and finding roots of the characteristic equation, which arises from a linearization of the model equations. Numerical simulations of the full model equations were conducted to assist in the interpretation of the bifurcation analysis. These techniques revealed a complex parameter region that allows a variety of qualitatively different model solutions: a regime having one stable, time-independent steady-state solution; regimes having one stable oscillatory solution only; and regimes having multiple possible stable oscillatory solutions. Model results suggest that the compliance of the TAL walls increases the tendency of the model TGF system to oscillate.     

Computational model of excitatory/inhibitory ratio imbalance role in attention deficit disorders

Impairments in attentional behaviors, including over-selectivity, under-selectivity, distractibility and difficulty in shift of attention, are widely reported in several developmental disorders, including autism. Uncharacteristic inhibitory to excitatory neuronal number ratio (IER) and abnormal synaptic strength levels in the brain are two broadly accepted neurobiological disorders observed in autistic individuals. These neurobiological findings are contrasting and their relation to the atypical attentional behaviors is not clear yet. In this paper, we take a computational approach to investigate the relation of imbalanced IER and abnormal synaptic strength to some well-documented spectrum of attentional impairments. The computational model is based on a modified version of a biologically plausible neural model of two competing minicolumns in IT cortex augmented with a simple model of top-down attention. Top-down attention is assumed to amplify (attenuates) attended (unattended) stimulus. The inhibitory synaptic strength parameter in the model is set such that typical attentional behavior is emerged. Then, according to related findings, the parameter is changed and the model's attentional behavior is considered. The simulation results show that, without any change in top-down attention, the abnormal inhibitory synaptic strength values - and IER imbalance- result in over-selectivity, under-selectivity, distractibility and difficulty in shift of attention in the model. It suggests that the modeled neurobiological abnormalities can be accounted for the attentional deficits. In addition, the atypical attentional behaviors do not necessarily point to impairments in top-down attention. Our simulations suggest that limited changes in the inhibitory synaptic strength and variations in top-down attention signal affect the model's attentional behaviors in the same way. So, limited deficits in the inhibitory strength may be alleviated by appropriate change in top-down attention biasing. Nevertheless, our model proposes that this compensation is not possible for very high and very low values of the inhibitory strength.     

Evolutionary optimization with data collocation for reverse engineering of biological networks

MOTIVATION: Modern experimental biology is moving away from analyses of single elements to whole-organism measurements. Such measured time-course data contain a wealth of information about the structure and dynamic of the pathway or network. The dynamic modeling of the whole systems is formulated as a reverse problem that requires a well-suited mathematical model and a very efficient computational method to identify the model structure and parameters. Numerical integration for differential equations and finding global parameter values are still two major challenges in this field of the parameter estimation of nonlinear dynamic biological systems. RESULTS: We compare three techniques of parameter estimation for nonlinear dynamic biological systems. In the proposed scheme, the modified collocation method is applied to convert the differential equations to the system of algebraic equations. The observed time-course data are then substituted into the algebraic system equations to decouple system interactions in order to obtain the approximate model profiles. Hybrid differential evolution (HDE) with population size of five is able to find a global solution. The method is not only suited for parameter estimation but also can be applied for structure identification. The solution obtained by HDE is then used as the starting point for a local search method to yield the refined estimates.     

Stability switches,oscillatory multistability,and spatio-temporal patterns of nonlinear oscillations in recurrently delay coupled neural networks

A model of time-delay recurrently coupled spatially segregated neural assemblies is here proposed. We show that it operates like some of the hierarchical architectures of the brain. Each assembly is a neural network with no delay in the local couplings between the units. The delay appears in the long range feedforward and feedback inter-assemblies communications. Bifurcation analysis of a simple four-units system in the autonomous case shows the richness of the dynamical behaviors in a biophysically plausible parameter region. We find oscillatory multistability, hysteresis, and stability switches of the rest state provoked by the time delay. Then we investigate the spatio-temporal patterns of bifurcating periodic solutions by using the symmetric local Hopf bifurcation theory of delay differential equations and derive the equation describing the flow on the center manifold that enables us determining the direction of Hopf bifurcations and stability of the bifurcating periodic orbits. We also discuss computational properties of the system due to the delay when an external drive of the network mimicks external sensory input.     

An accurate two-phase approximate solution to an acute viral infection model

During an acute viral infection, virus levels rise, reach a peak and then decline. Data and numerical solutions suggest the growth and decay phases are linear on a log scale. While viral dynamic models are typically nonlinear with analytical solutions difficult to obtain, the exponential nature of the solutions suggests approximations can be found. We derive a two-phase approximate solution to the target cell limited influenza model and illustrate its accuracy using data and previously established parameter values of six patients infected with influenza A. For one patient, the fall in virus concentration from its peak was not consistent with our predictions during the decay phase and an alternate approximation is derived. We find expressions for the rate and length of initial viral growth in terms of model parameters, the extent each parameter is involved in viral peaks, and the single parameter responsible for virus decay. We discuss applications of this analysis in antiviral treatments and in investigating host and virus heterogeneities.     

RankAggreg, an R package for weighted rank aggregation

Background

Determining the parameters of a mathematical model from quantitative measurements is the main bottleneck of modelling biological systems. Parameter values can be estimated from steady-state data or from dynamic data. The nature of suitable data for these two types of estimation is rather different. For instance, estimations of parameter values in pathway models, such as kinetic orders, rate constants, flux control coefficients or elasticities, from steady-state data are generally based on experiments that measure how a biochemical system responds to small perturbations around the steady state. In contrast, parameter estimation from dynamic data requires time series measurements for all dependent variables. Almost no literature has so far discussed the combined use of both steady-state and transient data for estimating parameter values of biochemical systems.

Results

In this study we introduce a constrained optimization method for estimating parameter values of biochemical pathway models using steady-state information and transient measurements. The constraints are derived from the flux connectivity relationships of the system at the steady state. Two case studies demonstrate the estimation results with and without flux connectivity constraints. The unconstrained optimal estimates from dynamic data may fit the experiments well, but they do not necessarily maintain the connectivity relationships. As a consequence, individual fluxes may be misrepresented, which may cause problems in later extrapolations. By contrast, the constrained estimation accounting for flux connectivity information reduces this misrepresentation and thereby yields improved model parameters.

Conclusion

The method combines transient metabolic profiles and steady-state information and leads to the formulation of an inverse parameter estimation task as a constrained optimization problem. Parameter estimation and model selection are simultaneously carried out on the constrained optimization problem and yield realistic model parameters that are more likely to hold up in extrapolations with the model.     

Dynamic model of hormonal systems coupled by negative feedback

Most hormone concentrations in the body are regulated by negative feedback mechanisms in which the production and release of hormones are regulated according to the concentration of related species. Also, it has been observed that several hormones are released in a variety of pulsatile patterns. In most cases, the mechanism driving these complex patterns is not well understood. Our model of two cells coupled through negative feedback to their external products demonstrates periodic, aperiodic and chaotic oscillations. The coupling between the cells seems to be responsible for these dynamic behaviors. The variety of dynamic behaviors observed in the model demonstrates that a simple physiological feedback loop mimicking the coupling between circulatory hormones and production centers could be the source of complex hormone release patterns observed in vivo.     

Diverse metabolic model parameters generate similar methionine cycle dynamics

Parameter estimation constitutes a major challenge in dynamic modeling of metabolic networks. Here we examine, via computational simulations, the influence of system nonlinearity and the nature of available data on the distribution and predictive capability of identified model parameters. Simulated methionine cycle metabolite concentration data (both with and without corresponding flux data) was inverted to identify model parameters consistent with it. Thousands of diverse parameter families were found to be consistent with the data to within moderate error, with most of the parameter values spanning over 1000-fold ranges irrespective of whether flux data was included. Due to strong correlations within the extracted parameter families, model predictions were generally reliable despite the broad ranges found for individual parameters. Inclusion of flux data, by strengthening these correlations, resulted in substantially more reliable flux predictions. These findings suggest that, despite the difficulty of extracting biochemically accurate model parameters from system level data, such data may nevertheless prove adequate for driving the development of predictive dynamic metabolic models.     

Modeling the interrelations between the calcium oscillations and ER membrane potential oscillations

A refined electrochemical model accounting for intracellular calcium oscillations and their interrelations with oscillations of the potential difference across the membrane of the endoplasmic reticulum (ER) or other intracellular calcium stores is established. The ATP dependent uptake of Ca2+ from the cytosol into the ER, the Ca2+ release from the ER through channels following a calcium-induced calcium release mechanism, and a potential-dependent Ca2+ leak flux out of the ER are included in the model and described by plausible rate laws. The binding of calcium to specific proteins such as calmodulin is taken into account. The quasi-electroneutrality condition allows us to express the transmembrane potential in terms of the concentrations of cytosolic calcium and free binding sites on proteins, which are the two independent variables of the model. We include monovalent ions in the model, because they make up a considerable portion in the balance of electroneutrality. As the permeability of the endoplasmic membrane for these ions is much higher than that for calcium ions, we assume the former to be in Nernst equilibrium. A stability analysis of the steady-state solutions (which are unique or multiple depending on parameter values) is carried out and the Hopf bifurcation leading from stable steady states to self-sustained oscillations is analysed with the help of appropriate mathematical techniques. The oscillations obtained by numerical integration exhibit the typical spike-like shape found in experiments and reasonable values of frequency and amplitude. The model describes the process of switching between stationary and pulsatile regimes as well as changes in oscillation frequency upon parameter changes. It turns out that calcium oscillations can arise without a permanent influx of calcium into the cell, when a calcium-buffering system such as calmodulin is included.     

Hopf bifurcations in multiple-parameter space of the Hodgkin-Huxley equations I. Global organization of bistable periodic solutions

 The Hodgkin-Huxley equations (HH) are parameterized by a number of parameters and shows a variety of qualitatively different behaviors depending on the parameter values. We explored the dynamics of the HH for a wide range of parameter values in the multiple-parameter space, that is, we examined the global structure of bifurcations of the HH. Results are summarized in various two-parameter bifurcation diagrams with I _ext (externally applied DC current) as the abscissa and one of the other parameters as the ordinate. In each diagram, the parameter plane was divided into several regions according to the qualitative behavior of the equations. In particular, we focused on periodic solutions emerging via Hopf bifurcations and identified parameter regions in which either two stable periodic solutions with different amplitudes and periods and a stable equilibrium point or two stable periodic solutions coexist. Global analysis of the bifurcation structure suggested that generation of these regions is associated with degenerate Hopf bifurcations. Received: 23 April 1999 / Accepted in revised form: 24 September 1999     

Improving evolutionary models of protein interaction networks

MOTIVATION: Theoretical models of biological networks are valuable tools in evolutionary inference. Theoretical models based on gene duplication and divergence provide biologically plausible evolutionary mechanics. Similarities found between empirical networks and their theoretically generated counterpart are considered evidence of the role modeled mechanics play in biological evolution. However, the method by which these models are parameterized can lead to questions about the validity of the inferences. Selecting parameter values in order to produce a particular topological value obfuscates the possibility that the model may produce a similar topology for a large range of parameter values. Alternately, a model may produce a large range of topologies, allowing (incorrect) parameter values to produce a valid topology from an otherwise flawed model. In order to lend biological credence to the modeled evolutionary mechanics, parameter values should be derived from the empirical data. Furthermore, recent work indicates that the timing and fate of gene duplications are critical to proper derivation of these parameters. RESULTS: We present a methodology for deriving evolutionary rates from empirical data that is used to parameterize duplication and divergence models of protein interaction network evolution. Our method avoids shortcomings of previous methods, which failed to consider the effect of subsequent duplications. From our parameter values, we find that concurrent and existing existing duplication and divergence models are insufficient for modeling protein interaction network evolution. We introduce a model enhancement based on heritable interaction sites on the surface of a protein and find that it more closely reflects the high clustering found in the empirical network.     

Probing protein colloidal behavior in membrane‐based separation processes using spectrofluorometric Rayleigh scattering data

One of the primary problems in membrane‐based protein separation is membrane fouling. In this study we explored the feasibility of employing Rayleigh light scattering data from fluorescence studies combined with chemometric techniques to determine whether a correlation could be established with membrane fouling phenomena. Membrane flux was measured in a dead‐end UF filtration system and the effect of protein solution properties on the flux decline was systematically investigated. A variety of proteins were used as a test case in this study. In parallel, the colloidal behavior of the protein solutions was assessed by employing multiwavelength Rayleigh scattering measurements. To assess the usefulness of Rayleigh scattering measurements for probing the colloidal behavior of proteins, a protein solution of β‐lactoglobulin was used as a base‐case scenario. The colloidal behavior of different β‐lactoglobulin solutions was inferred based on published data for this protein, under identical solution conditions, where techniques other than Rayleigh scattering had been used. Using this approach, good agreement was observed between scattering data and the colloidal behavior of this protein. To test the hypothesis that a high degree of aggregation will lead to increased membrane fouling, filtration data was used to find whether the Rayleigh scattering intensity correlated with permeate flux changes. It was found that for protein solutions which were stable and did not aggregate, fouling was reduced and these solutions exhibited reduced Rayleigh scattering. When the aggregation behavior of the solution was favored, significant flux declines occurred and were highly correlated with increased Rayleigh scattering. © 2010 American Institute of Chemical Engineers Biotechnol. Prog., 2010     

Singular Value Decomposition of the Radial Distribution Function for Hard Sphere and Square Well Potentials

We compute the singular value decomposition of the radial distribution function for hard sphere, and square well solutions. We find that decomposes into a small set of basis vectors allowing for an extremely accurate representation at all interpolated densities and potential strengths. In addition, we find that the coefficient vectors describing the magnitude of each basis vector are well described by a low-order polynomial. We provide a program to calculate in this compact representation for the investigated parameter range.     

A parity-structured matrix model for tsetse populations

A matrix model is used to describe the dynamics of a population of female tsetse flies structured by parity (i.e., by the number of larvae laid). For typical parameter values, the intrinsic growth rate of the population is zero when the adult daily survival rate is 0.970, corresponding to an adult life expectancy of 1/0.030 = 33.3 days. This value is plausible and consistent with results found earlier by others. The intrinsic growth rate is insensitive to the variance of the interlarval period. Temperature being a function of the time of the year, a known relationship between temperature and mean pupal and interlarval times was used to produce a time-varying version of the model which was fitted to temperature and (estimated) population data. With well-chosen parameter values, the modeled population replicated at least roughly the population data. This illustrates dynamically the abiotic effect of temperature on population growth. Given that tsetse flies are the vectors of trypanosomiasis ("sleeping sickness") the model provides a framework within which future transmission models can be developed in order to study the impact of altered temperatures on the spread of this deadly disease.     

Integrating behavioral and neural data in a model of zebrafish network interaction

The spinal neural networks of larval zebrafish (Danio rerio) generate a variety of movements such as escape, struggling, and swimming. Various mechanisms at the neural and network levels have been proposed to account for switches between these behaviors. However, there are currently no detailed demonstrations of such mechanisms. This makes determining which mechanisms are plausible extremely difficult. In this paper, we propose a detailed biologically plausible model of the interactions between the swimming and escape networks in the larval zebrafish, while taking into account anatomical and physiological evidence. We show that the results of our neural model generate the expected behavior when used to control a hydrodynamic model of carangiform locomotion. As a result, the model presented here is a clear demonstration of a plausible mechanism by which these distinct behaviors can be controlled. Interestingly, the networks are anatomically overlapping, despite clear differences in behavioral function and physiology.     

Robustness,flexibility, and the role of lateral inhibition in the neurogenic network

BACKGROUND: Many gene networks used by developing organisms have been conserved over long periods of evolutionary time. Why is that? We showed previously that a model of the segment polarity network in Drosophila is robust to parameter variation and is likely to act as a semiautonomous patterning module. Is this true of other networks as well? RESULTS: We present a model of the core neurogenic network in Drosophila. Our model exhibits at least three related pattern-resolving behaviors that the real neurogenic network accomplishes during embryogenesis in Drosophila. Furthermore, we find that it exhibits these behaviors across a wide range of parameter values, with most of its parameters able to vary more than an order of magnitude while it still successfully forms our test patterns. With a single set of parameters, different initial conditions (prepatterns) can select between different behaviors in the network's repertoire. We introduce two new measures for quantifying network robustness that mimic recombination and allelic divergence and use these to reveal the shape of the domain in the parameter space in which the model functions. We show that lateral inhibition yields robustness to changes in prepatterns and suggest a reconciliation of two divergent sets of experimental results. Finally, we show that, for this model, robustness confers functional flexibility. CONCLUSIONS: The neurogenic network is robust to changes in parameter values, which gives it the flexibility to make new patterns. Our model also offers a possible resolution of a debate on the role of lateral inhibition in cell fate specification.