Characterizing multistationarity regimes in biochemical reaction networks

Switch like responses appear as common strategies in the regulation of cellular systems. Here we present a method to characterize bistable regimes in biochemical reaction networks that can be of use to both direct and reverse engineering of biological switches. In the design of a synthetic biological switch, it is important to study the capability for bistability of the underlying biochemical network structure. Chemical Reaction Network Theory (CRNT) may help at this level to decide whether a given network has the capacity for multiple positive equilibria, based on their structural properties. However, in order to build a working switch, we also need to ensure that the bistability property is robust, by studying the conditions leading to the existence of two different steady states. In the reverse engineering of biological switches, knowledge collected about the bistable regimes of the underlying potential model structures can contribute at the model identification stage to a drastic reduction of the feasible region in the parameter space of search. In this work, we make use and extend previous results of the CRNT, aiming not only to discriminate whether a biochemical reaction network can exhibit multiple steady states, but also to determine the regions within the whole space of parameters capable of producing multistationarity. To that purpose we present and justify a condition on the parameters of biochemical networks for the appearance of multistationarity, and propose an efficient and reliable computational method to check its satisfaction through the parameter space.     

Robustness analysis and behavior discrimination in enzymatic reaction networks

Characterizing the behavior and robustness of enzymatic networks with numerous variables and unknown parameter values is a major challenge in biology, especially when some enzymes have counter-intuitive properties or switch-like behavior between activation and inhibition. In this paper, we propose new methodological and tool-supported contributions, based on the intuitive formalism of temporal logic, to express in a rigorous manner arbitrarily complex dynamical properties. Our multi-step analysis allows efficient sampling of the parameter space in order to define feasible regions in which the model exhibits imposed or experimentally observed behaviors. In a first step, an algorithmic methodology involving sensitivity analysis is conducted to determine bifurcation thresholds for a limited number of model parameters or initial conditions. In a second step, this boundary detection is supplemented by a global robustness analysis, based on quasi-Monte Carlo approach that takes into account all model parameters. We apply this method to a well-documented enzymatic reaction network describing collagen proteolysis by matrix metalloproteinase MMP2 and membrane type 1 metalloproteinase (MT1-MMP) in the presence of tissue inhibitor of metalloproteinase TIMP2. For this model, our method provides an extended analysis and quantification of network robustness toward paradoxical TIMP2 switching activity between activation or inhibition of MMP2 production. Further implication of our approach is illustrated by demonstrating and analyzing the possible existence of oscillatory behaviors when considering an extended open configuration of the enzymatic network. Notably, we construct bifurcation diagrams that specify key parameters values controlling the co-existence of stable steady and non-steady oscillatory proteolytic dynamics.     

Case study and application of process analytical technology (PAT) towards bioprocessing: Use of tryptophan fluorescence as at‐line tool for making pooling decisions for process chromatography

Process analytical technology (PAT) has been gaining momentum in the biopharmaceutical community due to the potential for continuous real time quality assurance resulting in improved operational control and compliance. Two imperatives for implementing any PAT tool are that “variability is managed by the process” and “product quality attributes can be accurately and reliably predicted over the design space established for materials used, process parameters, manufacturing, environmental, and other conditions.” Recently, we have been examining the feasibility of applying different analytical tools to bioprocessing unit operations. We have previously demonstarted that commercially available online‐high performance liquid chromatography and ultra performance liquid chromatography systems can be used for analysis that can facilitate real‐time decisions for column pooling based on product quality attributes (Rathore et al., 2008 a,b). In this article, we review an at‐line tool that can be used for pooling of process chromatography columns. We have demonstrated that our tryptophan fluorescence method offers a feasible approach and meets the requirements of a PAT application. It is significantly faster than the alternative of fractionation, offline analysis followed by pooling. Although the method as presented here is not an online method, this technique may offer better resolution for certain applications and may be a more optimal approach as it is very conducive to implementation in a manufacturing environment. This technique is also amenable to be used as an online tool via front face fluorescence measurements done concurrently with product concentration determination by UV. © 2009 American Institute of Chemical Engineers Biotechnol. Prog., 2009     

Allostery and conformational free energy changes in human tryptophanyl‐tRNA synthetase from essential dynamics and structure networks

The interdependence of the concept of allostery and enzymatic catalysis, and they being guided by conformational mobility is gaining increased prominence. However, to gain a molecular level understanding of allostery and hence of enzymatic catalysis, it is of utter importance that the networks of amino acids participating in allostery be deciphered. Our lab has been exploring the methods of network analysis combined with molecular dynamics simulations to understand allostery at molecular level. Earlier we had outlined methods to obtain communication paths and then to map the rigid/flexible regions of proteins through network parameters like the shortest correlated paths, cliques, and communities. In this article, we advance the methodology to estimate the conformational populations in terms of cliques/communities formed by interactions including the side‐chains and then to compute the ligand‐induced population shift. Finally, we obtain the free‐energy landscape of the protein in equilibrium, characterizing the free‐energy minima accessed by the protein complexes. We have chosen human tryptophanyl‐tRNA synthetase (hTrpRS), a protein responsible for charging tryptophan to its cognate tRNA during protein biosynthesis for this investigation. This is a multidomain protein exhibiting excellent allosteric communication. Our approach has provided valuable structural as well as functional insights into the protein. The methodology adopted here is highly generalized to illuminate the linkage between protein structure networks and conformational mobility involved in the allosteric mechanism in any protein with known structure. Proteins 2010. © 2009 Wiley‐Liss, Inc.     

From Boolean Network Model to Continuous Model Helps in Design of Functional Circuits

Computational circuit design with desired functions in a living cell is a challenging task in synthetic biology. To achieve this task, numerous methods that either focus on small scale networks or use evolutionary algorithms have been developed. Here, we propose a two-step approach to facilitate the design of functional circuits. In the first step, the search space of possible topologies for target functions is reduced by reverse engineering using a Boolean network model. In the second step, continuous simulation is applied to evaluate the performance of these topologies. We demonstrate the usefulness of this method by designing an example biological function: the SOS response of E. coli. Our numerical results show that the desired function can be faithfully reproduced by candidate networks with different parameters and initial conditions. Possible circuits are ranked according to their robustness against perturbations in parameter and gene expressions. The biological network is among the candidate networks, yet novel designs can be generated. Our method provides a scalable way to design robust circuits that can achieve complex functions, and makes it possible to uncover design principles of biological networks.     

Model‐based identifiable parameter determination applied to a simultaneous saccharification and fermentation process model for bio‐ethanol production

In this work, a methodology for the model‐based identifiable parameter determination (MBIPD) is presented. This systematic approach is proposed to be used for structure and parameter identification of nonlinear models of biological reaction networks. Usually, this kind of problems are over‐parameterized with large correlations between parameters. Hence, the related inverse problems for parameter determination and analysis are mathematically ill‐posed and numerically difficult to solve. The proposed MBIPD methodology comprises several tasks: (i) model selection, (ii) tracking of an adequate initial guess, and (iii) an iterative parameter estimation step which includes an identifiable parameter subset selection (SsS) algorithm and accuracy analysis of the estimated parameters. The SsS algorithm is based on the analysis of the sensitivity matrix by rank revealing factorization methods. Using this, a reduction of the parameter search space to a reasonable subset, which can be reliably and efficiently estimated from available measurements, is achieved. The simultaneous saccharification and fermentation (SSF) process for bio‐ethanol production from cellulosic material is used as case study for testing the methodology. The successful application of MBIPD to the SSF process demonstrates a relatively large reduction in the identified parameter space. It is shown by a cross‐validation that using the identified parameters (even though the reduction of the search space), the model is still able to predict the experimental data properly. Moreover, it is shown that the model is easily and efficiently adapted to new process conditions by solving reduced and well conditioned problems. © 2013 American Institute of Chemical Engineers Biotechnol. Prog., 29:1064–1082, 2013     

Robust parameter estimation during logistic modeling of batch and fed‐batch culture kinetics

Methods for robust logistic modeling of batch and fed‐batch mammalian cell cultures are presented in this study. Linearized forms of the logistic growth, logistic decline, and generalized logistic equation were derived to obtain initial estimates of the parameters by linear least squares. These initial estimates facilitated subsequent determination of refined values by nonlinear optimization using three different algorithms. Data from BHK, CHO, and hybridoma cells in batch or fed‐batch cultures at volumes ranging from 100 mL–300 L were tested with the above approach and solution convergence was obtained for all three nonlinear optimization approaches for all data sets. This result, despite the sensitivity of logistic equations to parameter variation because of their exponential nature, demonstrated that robust estimation of logistic parameters was possible by this combination of linearization followed by nonlinear optimization. The approach is relatively simple and can be implemented in a spreadsheet to robustly model mammalian cell culture batch or fed‐batch data. © 2009 American Institute of Chemical Engineers Biotechnol. Prog., 2009     

How affinity influences tolerance in an idiotypic network

Idiotypic network models give one possible justification for the appearance of tolerance for a certain category of cells while maintaining immunization for the others. In this paper, we provide new evidence that the manner in which affinity is defined in an idiotypic network model imposes a definite topology on the connectivity of the potential idiotypic network that can emerge. The resulting topology is responsible for very different qualitative behaviour of the network. We show that using a 2D shape-space model with affinity based on complementary regions, a cluster-free topology results that clearly divides the space into distinct zones; if antigens fall into a zone in which there are no available antibodies to bind to, they are tolerated. On the other hand, if they fall into a zone in which there are highly concentrated antibodies available for binding, then they will be eliminated. On the contrary, using a 2D shape space with an affinity function based on cell similarity, a highly clustered topology emerges in which there is no separation of the space into isolated tolerant and non-tolerant zones. Using a bit-string shape space, both similar and complementary affinity measures also result in highly clustered networks. In the networks whose topologies exhibit high clustering, the tolerant and intolerant zones are so intertwined that the networks either reject all antigen or tolerate all antigen. We show that the distribution and topology of the antibody network defined by the complete set of nodes and links-an autonomous feature of the system-therefore selects which antigens are tolerated and which are eliminated.     

Geometry and topology of parameter space: investigating measures of robustness in regulatory networks

The concept of robustness of regulatory networks has been closely related to the nature of the interactions among genes, and the capability of pattern maintenance or reproducibility. Defining this robustness property is a challenging task, but mathematical models have often associated it to the volume of the space of admissible parameters. Not only the volume of the space but also its topology and geometry contain information on essential aspects of the network, including feasible pathways, switching between two parallel pathways or distinct/disconnected active regions of parameters. A method is presented here to characterize the space of admissible parameters, by writing it as a semi-algebraic set, and then theoretically analyzing its topology and geometry, as well as volume. This method provides a more objective and complete measure of the robustness of a developmental module. As a detailed case study, the segment polarity gene network is analyzed.     

Elementary Flux Modes Analysis of Functional Domain Networks Allows a Better Metabolic Pathway Interpretation

Metabolic network analysis is an important step for the functional understanding of biological systems. In these networks, enzymes are made of one or more functional domains often involved in different catalytic activities. Elementary flux mode (EFM) analysis is a method of choice for the topological studies of these enzymatic networks. In this article, we propose to use an EFM approach on networks that encompass available knowledge on structure-function. We introduce a new method that allows to represent the metabolic networks as functional domain networks and provides an application of the algorithm for computing elementary flux modes to analyse them. Any EFM that can be represented using the classical representation can be represented using our functional domain network representation but the fine-grained feature of functional domain networks allows to highlight new connections in EFMs. This methodology is applied to the tricarboxylic acid cycle (TCA cycle) of Bacillus subtilis, and compared to the classical analyses. This new method of analysis of the functional domain network reveals that a specific inhibition on the second domain of the lipoamide dehydrogenase (pdhD) component of pyruvate dehydrogenase complex leads to the loss of all fluxes. Such conclusion was not predictable in the classical approach.     

Dose-response aligned circuits in signaling systems

Cells use biological signal transduction pathways to respond to environmental stimuli and the behavior of many cell types depends on precise sensing and transmission of external information. A notable property of signal transduction that was characterized in the Saccharomyces cerevisiae yeast cell and many mammalian cells is the alignment of dose-response curves. It was found that the dose response of the receptor matches closely the dose responses of the downstream. This dose-response alignment (DoRA) renders equal sensitivities and concordant responses in different parts of signaling system and guarantees a faithful information transmission. The experimental observations raise interesting questions about the nature of the information transmission through DoRA signaling networks and design principles of signaling systems with this function. Here, we performed an exhaustive computational analysis on network architectures that underlie the DoRA function in simple regulatory networks composed of two and three enzymes. The minimal circuits capable of DoRA were examined with Michaelis-Menten kinetics. Several motifs that are essential for the dynamical function of DoRA were identified. Systematic analysis of the topology space of robust DoRA circuits revealed that, rather than fine-tuning the network's parameters, the function is primarily realized by enzymatic regulations on the controlled node that are constrained in limiting regions of saturation or linearity.     

Automatic network fingerprinting through single-node motifs

Complex networks have been characterised by their specific connectivity patterns (network motifs), but their building blocks can also be identified and described by node-motifs-a combination of local network features. One technique to identify single node-motifs has been presented by Costa et al. (L. D. F. Costa, F. A. Rodrigues, C. C. Hilgetag, and M. Kaiser, Europhys. Lett., 87, 1, 2009). Here, we first suggest improvements to the method including how its parameters can be determined automatically. Such automatic routines make high-throughput studies of many networks feasible. Second, the new routines are validated in different network-series. Third, we provide an example of how the method can be used to analyse network time-series. In conclusion, we provide a robust method for systematically discovering and classifying characteristic nodes of a network. In contrast to classical motif analysis, our approach can identify individual components (here: nodes) that are specific to a network. Such special nodes, as hubs before, might be found to play critical roles in real-world networks.     

Noise propagation in two-step series MAPK cascade

Series MAPK enzymatic cascades, ubiquitously found in signaling networks, act as signal amplifiers and play a key role in processing information during signal transduction in cells. In activated cascades, cell-to-cell variability or noise is bound to occur and thereby strongly affects the cellular response. Commonly used linearization method (LM) applied to Langevin type stochastic model of the MAPK cascade fails to accurately predict intrinsic noise propagation in the cascade. We prove this by using extensive stochastic simulations for various ranges of biochemical parameters. This failure is due to the fact that the LM ignores the nonlinear effects on the noise. However, LM provides a good estimate of the extrinsic noise propagation. We show that the correct estimate of intrinsic noise propagation in signaling networks that contain at least one enzymatic step can be obtained only through stochastic simulations. Noise propagation in the cascade depends on the underlying biochemical parameters which are often unavailable. Based on a combination of global sensitivity analysis (GSA) and stochastic simulations, we developed a systematic methodology to characterize noise propagation in the cascade. GSA predicts that noise propagation in MAPK cascade is sensitive to the total number of upstream enzyme molecules and the total number of molecules of the two substrates involved in the cascade. We argue that the general systematic approach proposed and demonstrated on MAPK cascade must accompany noise propagation studies in biological networks.     

Computational biotechnology: Prediction of competitive substrate inhibition of enzymes by buffer compounds with protein-ligand docking

In vitro enzymatic activity highly depends on the reaction medium. One of the most important parameters is the buffer used to keep the pH stable. The buffering compound prevents a severe pH-change and therefore a possible denaturation of the enzyme. However buffer agents can also have negative effects on the enzymatic activity, such as competitive substrate inhibition. We assess this effect with a computational approach based on a protein-ligand docking method and the HYDE scoring function. Our method predicts competitive binding of the buffer compound to the active site of the enzyme. Using data from literature and new experimental data, the procedure is evaluated on nine different enzymatic reactions. The method predicts buffer-enzyme interactions and is able to score these interactions with the correct trend of enzymatic activities. Using the new method, possible buffers can be selected or discarded prior to laboratory experiments.     

Enzymatic grafting of hexyloxyphenol onto chitosan to alter surface and rheological properties

An enzymatic method to graft hexyloxyphenol onto the biopolymer chitosan was studied. The method employs tyrosinase to convert the phenol into a reactive o-quinone, which undergoes subsequent nonenzymatic reaction with chitosan. Reactions were conducted under heterogeneous conditions using chitosan films and also under homogeneous conditions using aqueous methanolic mixtures capable of dissolving both hexyloxyphenol and chitosan. Tyrosinase was shown to catalyze the oxidation of hexyloxyphenol in such aqueous methanolic solutions. Chemical evidence for covalent grafting onto chitosan was provided by three independent spectroscopic approaches. Specifically, enzymatic modification resulted in (1) the appearance of broad absorbance in the 350-nm region of the UV/vis spectra for chitosan films; (2) changes in the NH bending and stretching regions of chitosan's IR spectra; and (3) a base-soluble material with (1)H-NMR signals characteristic of both chitosan and the alkyl groups of hexyloxyphenol. Hexyloxyphenol modification resulted in dramatic changes in chitosan's functional properties. On the basis of contact angle measurements, heterogeneous modification of a chitosan film yielded a hydrophobic surface. Homogeneously modified chitosan offered rheological properties characteristic of associating water-soluble polymers.     

Enzyme immobilization via silaffin‐mediated autoencapsulation in a biosilica support

Enzymes and other biomolecules are often immobilized in a matrix to improve their stability or to improve their ability to be reused. Performing a polycondensation reaction in the presence of a biomolecule of interest relies on random entrapment events during polymerization and may not ensure efficient, homogeneous, or complete biomolecule encapsulation. To overcome these limitations, we have developed a method of incorporating autosilification activity into proteins without affecting enzymatic functionality. The unmodified R5 silaffin peptide from Cylindrotheca fusiformis is capable of initiating silica polycondensation in vitro at ambient temperatures and pressures in aqueous solution. In this study, translational fusion proteins between R5 and various functional proteins (phosphodiesterase, organophosphate hydrolase, and green fluorescent protein) were produced in Escherichia coli. Each of the fusion proteins initiated silica polycondensation, and enzymatic activity (or fluorescence) was retained in the resulting silica spheres. Under certain circumstances, the enzymatically‐active biosilica displayed improved stability relative to free enzyme at elevated temperatures. © 2009 American Institute of Chemical Engineers Biotechnol. Prog., 2009     

Exact Probability Distributions of Selected Species in Stochastic Chemical Reaction Networks

Chemical reactions are discrete, stochastic events. As such, the species’ molecular numbers can be described by an associated master equation. However, handling such an equation may become difficult due to the large size of reaction networks. A commonly used approach to forecast the behaviour of reaction networks is to perform computational simulations of such systems and analyse their outcome statistically. This approach, however, might require high computational costs to provide accurate results. In this paper we opt for an analytical approach to obtain the time-dependent solution of the Chemical Master Equation for selected species in a general reaction network. When the reaction networks are composed exclusively of zeroth and first-order reactions, this analytical approach significantly alleviates the computational burden required by simulation-based methods. By building upon these analytical solutions, we analyse a general monomolecular reaction network with an arbitrary number of species to obtain the exact marginal probability distribution for selected species. Additionally, we study two particular topologies of monomolecular reaction networks, namely (i) an unbranched chain of monomolecular reactions with and without synthesis and degradation reactions and (ii) a circular chain of monomolecular reactions. We illustrate our methodology and alternative ways to use it for non-linear systems by analysing a protein autoactivation mechanism. Later, we compare the computational load required for the implementation of our results and a pure computational approach to analyse an unbranched chain of monomolecular reactions. Finally, we study calcium ions gates in the sarco/endoplasmic reticulum mediated by ryanodine receptors.     

Shape, Size, and Robustness: Feasible Regions in the Parameter Space of Biochemical Networks

The concept of robustness of regulatory networks has received much attention in the last decade. One measure of robustness has been associated with the volume of the feasible region, namely, the region in the parameter space in which the system is functional. In this paper, we show that, in addition to volume, the geometry of this region has important consequences for the robustness and the fragility of a network. We develop an approximation within which we could algebraically specify the feasible region. We analyze the segment polarity gene network to illustrate our approach. The study of random walks in the parameter space and how they exit the feasible region provide us with a rich perspective on the different modes of failure of this network model. In particular, we found that, between two alternative ways of activating Wingless, one is more robust than the other. Our method provides a more complete measure of robustness to parameter variation. As a general modeling strategy, our approach is an interesting alternative to Boolean representation of biochemical networks.     

Evolutionary innovations and the organization of protein functions in genotype space

The organization of protein structures in protein genotype space is well studied. The same does not hold for protein functions, whose organization is important to understand how novel protein functions can arise through blind evolutionary searches of sequence space. In systems other than proteins, two organizational features of genotype space facilitate phenotypic innovation. The first is that genotypes with the same phenotype form vast and connected genotype networks. The second is that different neighborhoods in this space contain different novel phenotypes. We here characterize the organization of enzymatic functions in protein genotype space, using a data set of more than 30,000 proteins with known structure and function. We show that different neighborhoods of genotype space contain proteins with very different functions. This property both facilitates evolutionary innovation through exploration of a genotype network, and it constrains the evolution of novel phenotypes. The phenotypic diversity of different neighborhoods is caused by the fact that some functions can be carried out by multiple structures. We show that the space of protein functions is not homogeneous, and different genotype neighborhoods tend to contain a different spectrum of functions, whose diversity increases with increasing distance of these neighborhoods in sequence space. Whether a protein with a given function can evolve specific new functions is thus determined by the protein's location in sequence space.     

Deterministic epidemiological models at the individual level

In many fields of science including population dynamics, the vast state spaces inhabited by all but the very simplest of systems can preclude a deterministic analysis. Here, a class of approximate deterministic models is introduced into the field of epidemiology that reduces this state space to one that is numerically feasible. However, these reduced state space master equations do not in general form a closed set. To resolve this, the equations are approximated using closure approximations. This process results in a method for constructing deterministic differential equation models with a potentially large scope of application including dynamic directed contact networks and heterogeneous systems using time dependent parameters. The method is exemplified in the case of an SIR (susceptible-infectious-removed) epidemiological model and is numerically evaluated on a range of networks from spatially local to random. In the context of epidemics propagated on contact networks, this work assists in clarifying the link between stochastic simulation and traditional population level deterministic models.