Search

Abstract Single-cell omics technologies can measure millions of cells for up to thousands of biomolecular features, enabling data-driven studies of complex biological networks. However, these high-throughput experimental techniques often cannot track individual cells over time, thus complicating the understanding of dynamics such as time trajectories of cell states. These “dynamical phenotypes” are key to understanding biological phenomena such as differentiation fates. We show by mathematical analysis that, in spite of high dimensionality and lack of individual cell traces, three time-points of single-cell omics data are theoretically necessary and sufficient to uniquely determine the network interaction matrix and associated dynamics. Moreover, we show through numerical simulations that an interaction matrix can be accurately determined with three or more time-points even in the presence of sampling and measurement noise typical of single-cell omics. Our results can guide the design of single-cell omics time-course experiments, and provide a tool for data-driven phase-space analysis.

Hong, Kim and colleagues develop a method for analytically deriving the stationary distributions of stochastic biochemical reaction networks using network structure transformation. They provide this method in a user-friendly computational package, CASTANET, available open-source.

Cell-fate networks are traditionally studied within the framework of gene regulatory networks. This paradigm considers only interactions of genes through expressed transcription factors and does not incorporate chromatin modification processes. This paper introduces a mathematical model that seamlessly combines gene regulatory networks and DNA methylation (DNAm), with the goal of quantitatively characterizing the contribution of epigenetic regulation to gene silencing. The ‘Basin of Attraction percentage’ is introduced as a metric to quantify gene silencing abilities. As a case study, a computational and theoretical analysis is carried out for a model of the pluripotent stem cell circuit as well as a simplified self-activating gene model. The results confirm that the methodology quantitatively captures the key role that DNAm plays in enhancing the stability of the silenced gene state.

Summary: A design for genetically encoded counters is proposed via repressor-based circuits. An N-bit counter reads sequences of input pulses and displays the total number of pulses, modulo 2N. The design is based on distributed computation with specialized cell types allocated to specific tasks. This allows scalability and bypasses constraints on the maximal number of circuit genes per cell due to toxicity or failures due to resource limitations. The design starts with a single-bit counter. The N-bit counter is then obtained by interconnecting (using diffusible chemicals) a set of N single-bit counters and connector modules. An optimization framework is used to determine appropriate gate parameters and to compute bounds on admissible pulse widths and relaxation (inter-pulse) times, as well as to guide the construction of novel gates. This work can be viewed as a step toward obtaining circuits that are capable of finite automaton computation in analogy to digital central processing units.

A dynamical system entrains to a periodic input if its state converges globally to an attractor with the same period. In particular, for a constant input, the state converges to a unique equilibrium point for any initial condition. We consider the problem of maximizing a weighted average of the system’s output along the periodic attractor. The gain of entrainment is the benefit achieved by using a non-constant periodic input relative to a constant input with the same time average. Such a problem amounts to optimal allocation of resources in a periodic manner. We formulate this problem as a periodic optimal control problem, which can be analysed by means of the Pontryagin maximum principle or solved numerically via powerful software packages. We then apply our framework to a class of nonlinear occupancy models that appear frequently in biological synthesis systems and other applications. We show that, perhaps surprisingly, constant inputs are optimal for various architectures. This suggests that the presence of non-constant periodic signals, which frequently appear in biological occupancy systems, is a signature of an underlying time-varying objective functional being optimized.

Exact analytical and closed-form solutions to the transient diffusion in bi-layer composites with external mass transfer resistance are reported. Expressions for the concentrations and the mass permeated are derived in both the Laplace and time domains through the use of the Laplace transform Inversion Theorem. The lead and lag times, which are often of importance in the characterization of membranes and arise from the analysis of the asymptotic behavior of the mass permeated through the bi-layer composite, were also derived. The presented solutions are also compared to previously derived limiting cases of the diffusion in a bi-layer with an impermeable wall and constant concentrations at the upstream and downstream boundaries. Analysis of the time lag shows that this membrane property is independent of the direction of flow. Finally, an outline is provided of how these transient solutions in response to a step function increase in concentration can be used to derive more complex input conditions. The importance of adequately handling boundary layer effects has a wide array of applications such as the study of bi-layers undergoing phenomena of heat convection, gas film resistance, and absorption/desorption.

MotivationCirculating tumor cells (CTCs) are widely studied using liquid biopsy methods that analyze fractionally-small peripheral blood (PB) samples. However, little is known about natural fluctuations in CTC numbers that may occur over short timescales in vivo, and how these may affect detection and enumeration of rare CTCs from small blood samples.MethodsWe recently developed an optical instrument called “diffuse in vivo flow cytometry” (DiFC) that uniquely allows continuous, non-invasive counting of rare, green fluorescent protein expressing CTCs in large blood vessels in mice. Here, we used DiFC to study short-term changes in CTC numbers in multiple myeloma and Lewis lung carcinoma xenograft models. We analyzed CTC detections in over 100 h of DiFC data, and considered intervals corresponding to approximately 1%, 5%, 10%, and 20% of the PB volume. In addition, we analyzed changes in CTC numbers over 24 h (diurnal) periods.ResultsFor rare CTCs (fewer than 1 CTC per ml of blood), the use of short DiFC intervals (corresponding to small PB samples) frequently resulted in no detections. For more abundant CTCs, CTC numbers frequently varied by an order of magnitude or more over the time-scales considered. This variance in CTC detections far exceeded that expected by Poisson statistics or by instrument variability. Rather, the data were consistent with significant changes in mean numbers of CTCs on the timescales of minutes and hours.ConclusionsThe observed temporal changes can be explained by known properties of CTCs, namely, the continuous shedding of CTCs from tumors and the short half-life of CTCs in blood. It follows that the number of cells in a blood sample are strongly impacted by the timing of the draw. The issue is likely to be compounded for multicellular CTC clusters or specific CTC subtypes, which are even more rare than single CTCs. However, we show that enumeration can in principle be improved by averaging multiple samples, analysis of larger volumes, or development of methods for enumeration of CTCs directly in vivo.

One of the most important factors limiting the success of chemotherapy in cancer treatment is the phenomenon of drug resistance. We have recently introduced a framework for quantifying the effects of induced and non-induced resistance to cancer chemotherapy (Greene et al., 2018a, 2019). In this work, we expound on the details relating to an optimal control problem outlined in Greene et al. (2018a). The control structure is precisely characterized as a concatenation of bang-bang and path-constrained arcs via the Pontryagin Maximum Principle and differential Lie algebraic techniques. A structural identifiability analysis is also presented, demonstrating that patient-specific parameters may be measured and thus utilized in the design of optimal therapies prior to the commencement of therapy. For completeness, a detailed analysis of existence results is also included.

Metronomic chemotherapy can drastically enhance immunogenic tumor cell death. However, the mechanisms responsible are still incompletely understood. Here, we develop a mathematical model to elucidate the underlying complex interactions between tumor growth, immune system activation, and therapy-mediated immunogenic cell death. Our model is conceptually simple, yet it provides a surprisingly excellent fit to empirical data obtained from a GL261 SCID mouse glioma model treated with cyclophosphamide on a metronomic schedule. The model includes terms representing immune recruitment as well as the emergence of drug resistance during prolonged metronomic treatments. Strikingly, a single fixed set of parameters, adjusted neither for individuals nor for drug schedule, recapitulates experimental data across various drug regimens remarkably well, including treatments administered at intervals ranging from 6 to 12 days. Additionally, the model predicts peak immune activation times, rediscovering experimental data that had not been used in parameter fitting or in model construction. Notably, the validated model suggests that immunostimulatory and immunosuppressive intermediates are responsible for the observed phenomena of resistance and immune cell recruitment, and thus for variation of responses with respect to different schedules of drug administration.

It was recently reported that acute influenza infection of the lung promoted distal melanoma growth in the dermis of mice. Melanoma-specific CD8+ T cells were shunted to the lung in the presence of the infection, where they expressed high levels of inflammation-induced cell-activation blocker PD-1, and became incapable of migrating back to the tumor site. At the same time, co-infection virus-specific CD8+ T cells remained functional while the infection was cleared. It was also unexpectedly found that PD-1 blockade immunotherapy reversed this effect. Here, we proceed to ground the experimental observations in a mechanistic immunobiochemical model that incorporates T cell pathways that control PD-1 expression. A core component of our model is a kinetic motif, which we call a PD-1 Double Incoherent Feed-Forward Loop (DIFFL), and which reflects known interactions between IRF4, Blimp-1, and Bcl-6. The different activity levels of the PD-1 DIFFL components, as a function of the cognate antigen levels and the given inflammation context, manifest themselves in phenotypically distinct outcomes. Collectively, the model allowed us to put forward a few working hypotheses as follows: (i) the melanoma-specific CD8+ T cells re-circulating with the blood flow enter the lung where they express high levels of inflammation-induced cell-activation blocker PD-1 in the presence of infection; (ii) when PD-1 receptors interact with abundant PD-L1, constitutively expressed in the lung, T cells loose motility; (iii) at the same time, virus-specific cells adapt to strong stimulation by their cognate antigen by lowering the transiently-elevated expression of PD-1, remaining functional and mobile in the inflamed lung, while the infection is cleared. The role that T cell receptor (TCR) activation and feedback loops play in the underlying processes are also highlighted and discussed. We hope that the results reported in our study could potentially contribute to the advancement of immunological approaches to cancer treatment and, as well, to a better understanding of a broader complexity of fundamental interactions between pathogens and tumors.

The aim of this paper is to study the qualitative behaviour predicted by a mathematical model for the initial stage of T-cell activation. The state variables in the model are the concentrations of phosphorylation states of the T-cell receptor (TCR) complex and the phosphatase SHP-1 in the cell. It is shown that these quantities cannot approach zero and that the model possesses more than one positive steady state for certain values of the parameters. It can also exhibit damped oscillations. It is proved that the chemical concentration which represents the degree of activation of the cell, that of the maximally phosphorylated form of the TCR complex, is, in general, a non-monotone function of the activating signal. In particular, there are cases where there is a value of the dissociation constant of the ligand from the receptor which produces a maximal activation of the T cell. This suggests that mechanisms taking place in the first few minutes after activation and included in the model studied in this paper suffice to explain the optimal dissociation time seen in experiments. In this way, the results of certain simulations in the literature have been confirmed rigorously and some important features which had not previously been seen have been discovered.

We study a single species in a chemostat, limited by two nutrients, and separate nutrient uptake from growth. For a broad class of uptake and growth functions it is proved that a nontrivial equilibrium may exist. Moreover, if it exists it is unique and globally stable, generalizing a result in [15]. [ABSTRACT FROM AUTHOR]

We prove the global asymptotic stability of a well-known delayed negative-feedback model of testosterone dynamics, which has been proposed as a model of oscillatory behavior. We establish stability (and hence the impossibility of oscillations) even in the presence of delays of arbitrary length. [ABSTRACT FROM AUTHOR]

The fundamental relaxation result for Lipschitz differential inclusions is the Filippov-Wazewski Relaxation Theorem, which provides approximations of trajectories of a relaxed inclusion on finite intervals. A complementary result is presented, which provides approximations on infinite intervals, but does not guarantee that the approximation and the reference trajectory satisfy the same initial condition. [ABSTRACT FROM AUTHOR]

Feedback mechanisms for synthetic gene circuits are necessary to provide robustness to external perturbations. Here the authors validate a biomolecular controller based on a sigma and anti-sigma factor to achieve stable gene expression in the face of external disturbances in an in vitro synthetic gene circuit.

Complex molecular biological processes such as transcription and translation, signal transduction, post-translational modification cascades, and metabolic pathways can be described in principle by biochemical reactions that explicitly take into account the sophisticated network of chemical interactions regulating cell life. The ability to deduce the possible qualitative behaviors of such networks from a set of reactions is a central objective and an ongoing challenge in the field of systems biology. Unfortunately, the construction of complete mathematical models is often hindered by a pervasive problem: despite the wealth of qualitative graphical knowledge about network interactions, the form of the governing nonlinearities and/or the values of kinetic constants are hard to uncover experimentally. The kinetics can also change with environmental variations. This work addresses the following question: given a set of reactions and without assuming a particular form for the kinetics, what can we say about the asymptotic behavior of the network? Specifically, it introduces a class of networks that are "structurally (mono) attractive" meaning that they are incapable of exhibiting multiple steady states, oscillation, or chaos by virtue of their reaction graphs. These networks are characterized by the existence of a universal energy-like function called a Robust Lyapunov function (RLF). To find such functions, a finite set of rank-one linear systems is introduced, which form the extremals of a linear convex cone. The problem is then reduced to that of finding a common Lyapunov function for this set of extremals. Based on this characterization, a computational package, Lyapunov-Enabled Analysis of Reaction Networks (LEARN), is provided that constructs such functions or rules out their existence. An extensive study of biochemical networks demonstrates that LEARN offers a new unified framework. Basic motifs, three-body binding, and genetic networks are studied first. The work then focuses on cellular signalling networks including various post-translational modification cascades, phosphotransfer and phosphorelay networks, T-cell kinetic proofreading, and ERK signalling. The Ribosome Flow Model is also studied.

The goal of many single-cell studies on eukaryotic cells is to gain insight into the biochemical reactions that control cell fate and state. In this paper we introduce the concept of Effective Stoichiometric Spaces (ESS) to guide the reconstruction of biochemical networks from multiplexed, fixed time-point, single-cell data. In contrast to methods based solely on statistical models of data, the ESS method leverages the power of the geometric theory of toric varieties to begin unraveling the structure of chemical reaction networks (CRN). This application of toric theory enables a data-driven mapping of covariance relationships in single-cell measurements into stoichiometric information, one in which each cell subpopulation has its associated ESS interpreted in terms of CRN theory. In the development of ESS we reframe certain aspects of the theory of CRN to better match data analysis. As an application of our approach we process cytomery- and image-based single-cell datasets and identify differences in cells treated with kinase inhibitors. Our approach is directly applicable to data acquired using readily accessible experimental methods such as Fluorescence Activated Cell Sorting (FACS) and multiplex immunofluorescence.

Phenotypical variability in the absence of genetic variation often reflects complex energetic landscapes associated with underlying gene regulatory networks (GRNs). In this view, different phenotypes are associated with alternative states of complex nonlinear systems: stable attractors in deterministic models or modes of stationary distributions in stochastic descriptions. We provide theoretical and practical characterizations of these landscapes, specifically focusing on stochastic Slow Promoter Kinetics (SPK), a time scale relevant when transcription factor binding and unbinding are affected by epigenetic processes like DNA methylation and chromatin remodeling. In this case, largely unexplored except for numerical simulations, adiabatic approximations of promoter kinetics are not appropriate. In contrast to the existing literature, we provide rigorous analytic characterizations of multiple modes. A general formal approach gives insight into the influence of parameters and the prediction of how changes in GRN wiring, for example through mutations or artificial interventions, impact the possible number, location, and likelihood of alternative states. We adapt tools from the mathematical field of singular perturbation theory to represent stationary distributions of Chemical Master Equations for GRNs as mixtures of Poisson distributions and obtain explicit formulas for the locations and probabilities of metastable states as a function of the parameters describing the system. As illustrations, the theory is used to tease out the role of cooperative binding in stochastic models in comparison to deterministic models, and applications are given to various model systems, such as toggle switches in isolation or in communicating populations, a synthetic oscillator, and a trans-differentiation network.

Biochemical reaction networks (BRNs) in a cell frequently consist of reactions with disparate timescales. The stochastic simulations of such multiscale BRNs are prohibitively slow due to high computational cost for the simulations of fast reactions. One way to resolve this problem uses the fact that fast species regulated by fast reactions quickly equilibrate to their stationary distribution while slow species are unlikely to be changed. Thus, on a slow timescale, fast species can be replaced by their quasi-steady state (QSS): their stationary conditional expectation values for given slow species. As the QSS are determined solely by the state of slow species, such replacement leads to a reduced model, where fast species are eliminated. However, it is challenging to derive the QSS in the presence of nonlinear reactions. While various approximation schemes for the QSS have been developed, they often lead to considerable errors. Here, we propose two classes of multiscale BRNs which can be reduced by deriving an exact QSS rather than approximations. Specifically, if fast species constitute either a feedforward network or a complex balanced network, the reduced model based on the exact QSS can be derived. Such BRNs are frequently observed in a cell as the feedforward network is one of fundamental motifs of gene or protein regulatory networks. Furthermore, complex balanced networks also include various types of fast reversible bindings such as bindings between transcriptional factors and gene regulatory sites. The reduced models based on exact QSS, which can be calculated by the computational packages provided in this work, accurately approximate the slow scale dynamics of the original full model with much lower computational cost.

A recent paper by Karin et al. introduced a mathematical notion called dynamical compensation (DC) of biological circuits. DC was shown to play an important role in glucose homeostasis as well as other key physiological regulatory mechanisms. Karin et al. went on to provide a sufficient condition to test whether a given system has the DC property. Here, we show how DC can be formulated in terms of a well-known concept in systems biology, statistics, and control theory-that of parameter structural non-identifiability. Viewing DC as a parameter identification problem enables one to take advantage of powerful theoretical and computational tools to test a system for DC. We obtain as a special case the sufficient criterion discussed by Karin et al. We also draw connections to system equivalence and to the fold-change detection property.

Abstract Synthetic genetic systems share resources with the host, including machinery for transcription and translation. Phage RNA polymerases (RNAPs) decouple transcription from the host and generate high expression. However, they can exhibit toxicity and lack accessory proteins (σ factors and activators) that enable switching between different promoters and modulation of activity. Here, we show that T7 RNAP (883 amino acids) can be divided into four fragments that have to be co‐expressed to function. The DNA‐binding loop is encoded in a C‐terminal 285‐aa ‘σ fragment’, and fragments with different specificity can direct the remaining 601‐aa ‘core fragment’ to different promoters. Using these parts, we have built a resource allocator that sets the core fragment concentration, which is then shared by multiple σ fragments. Adjusting the concentration of the core fragment sets the maximum transcriptional capacity available to a synthetic system. Further, positive and negative regulation is implemented using a 67‐aa N‐terminal ‘α fragment’ and a null (inactivated) σ fragment, respectively. The α fragment can be fused to recombinant proteins to make promoters responsive to their levels. These parts provide a toolbox to allocate transcriptional resources via different schemes, which we demonstrate by building a system which adjusts promoter activity to compensate for the difference in copy number of two plasmids.

Synthetic constructs in biotechnology, biocomputing, and modern gene therapy interventions are often based on plasmids or transfected circuits which implement some form of "on-off" switch. For example, the expression of a protein used for therapeutic purposes might be triggered by the recognition of a specific combination of inducers (e.g., antigens), and memory of this event should be maintained across a cell population until a specific stimulus commands a coordinated shut-off. The robustness of such a design is hampered by molecular ("intrinsic") or environmental ("extrinsic") noise, which may lead to spontaneous changes of state in a subset of the population and is reflected in the bimodality of protein expression, as measured for example using flow cytometry. In this context, a "majority-vote" correction circuit, which brings deviant cells back into the required state, is highly desirable, and quorum-sensing has been suggested as a way for cells to broadcast their states to the population as a whole so as to facilitate consensus. In this paper, we propose what we believe is the first such a design that has mathematically guaranteed properties of stability and auto-correction under certain conditions. Our approach is guided by concepts and theory from the field of "monotone" dynamical systems developed by M. Hirsch, H. Smith, and others. We benchmark our design by comparing it to an existing design which has been the subject of experimental and theoretical studies, illustrating its superiority in stability and self-correction of synchronization errors. Our stability analysis, based on dynamical systems theory, guarantees global convergence to steady states, ruling out unpredictable ("chaotic") behaviors and even sustained oscillations in the limit of convergence. These results are valid no matter what are the values of parameters, and are based only on the wiring diagram. The theory is complemented by extensive computational bifurcation analysis, performed for a biochemically-detailed and biologically-relevant model that we developed. Another novel feature of our approach is that our theorems on exponential stability of steady states for homogeneous or mixed populations are valid independently of the number N of cells in the population, which is usually very large (N ≫ 1) and unknown. We prove that the exponential stability depends on relative proportions of each type of state only. While monotone systems theory has been used previously for systems biology analysis, the current work illustrates its power for synthetic biology design, and thus has wider significance well beyond the application to the important problem of coordination of toggle switches.

Abstract Modularity plays a fundamental role in the prediction of the behavior of a system from the behavior of its components, guaranteeing that the properties of individual components do not change upon interconnection. Just as electrical, hydraulic, and other physical systems often do not display modularity, nor do many biochemical systems, and specifically, genetic networks. Here, we study the effect of interconnections on the input–output dynamic characteristics of transcriptional components, focusing on a property, which we call ‘retroactivity’, that plays a role analogous to non‐zero output impedance in electrical systems. In transcriptional networks, retroactivity is large when the amount of transcription factor is comparable to, or smaller than, the amount of promoter‐binding sites, or when the affinity of such binding sites is high. To attenuate the effect of retroactivity, we propose a feedback mechanism inspired by the design of amplifiers in electronics. We introduce, in particular, a mechanism based on a phosphorylation–dephosphorylation cycle. This mechanism enjoys a remarkable insulation property, due to the fast timescales of the phosphorylation and dephosphorylation reactions.

Periodic oscillations play an important role in many biomedical systems. Proper functioning of biological systems that respond to periodic signals requires the ability to synchronize with the periodic excitation. For example, the sleep/wake cycle is a manifestation of an internal timing system that synchronizes to the solar day. In the terminology of systems theory, the biological system must entrain or phase-lock to the periodic excitation. Entrainment is also important in synthetic biology. For example, connecting several artificial biological systems that entrain to a common clock may lead to a well-functioning modular system. The cell-cycle is a periodic program that regulates DNA synthesis and cell division. Recent biological studies suggest that cell-cycle related genes entrain to this periodic program at the gene translation level, leading to periodically-varying protein levels of these genes. The ribosome flow model (RFM) is a deterministic model obtained via a mean-field approximation of a stochastic model from statistical physics that has been used to model numerous processes including ribosome flow along the mRNA. Here we analyze the RFM under the assumption that the initiation and/or transition rates vary periodically with a common period T. We show that the ribosome distribution profile in the RFM entrains to this periodic excitation. In particular, the protein synthesis pattern converges to a unique periodic solution with period T. To the best of our knowledge, this is the first proof of entrainment in a mathematical model for translation that encapsulates aspects such as initiation and termination rates, ribosomal movement and interactions, and non-homogeneous elongation speeds along the mRNA. Our results support the conjecture that periodic oscillations in tRNA levels and other factors related to the translation process can induce periodic oscillations in protein levels, and may suggest a new approach for re-engineering genetic systems to obtain a desired, periodic, protein synthesis rate.

An ubiquitous property of biological sensory systems is adaptation: a step increase in stimulus triggers an initial change in a biochemical or physiological response, followed by a more gradual relaxation toward a basal, pre-stimulus level. Adaptation helps maintain essential variables within acceptable bounds and allows organisms to readjust themselves to an optimum and non-saturating sensitivity range when faced with a prolonged change in their environment. Recently, it was shown theoretically and experimentally that many adapting systems, both at the organism and single-cell level, enjoy a remarkable additional feature: scale invariance, meaning that the initial, transient behavior remains (approximately) the same even when the background signal level is scaled. In this work, we set out to investigate under what conditions a broadly used model of biochemical enzymatic networks will exhibit scale-invariant behavior. An exhaustive computational study led us to discover a new property of surprising simplicity and generality, uniform linearizations with fast output (ULFO), whose validity we show is both necessary and sufficient for scale invariance of three-node enzymatic networks (and sufficient for any number of nodes). Based on this study, we go on to develop a mathematical explanation of how ULFO results in scale invariance. Our work provides a surprisingly consistent, simple, and general framework for understanding this phenomenon, and results in concrete experimental predictions.

This paper addresses the problem of providing mathematical conditions that allow one to ensure that biological networks, such as transcriptional systems, can be globally entrained to external periodic inputs. Despite appearing obvious at first, this is by no means a generic property of nonlinear dynamical systems. Through the use of contraction theory, a powerful tool from dynamical systems theory, it is shown that certain systems driven by external periodic signals have the property that all their solutions converge to a fixed limit cycle. General results are proved, and the properties are verified in the specific cases of models of transcriptional systems as well as constructs of interest in synthetic biology. A self-contained exposition of all needed results is given in the paper.

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.

It has previously been shown that generic cortical microcircuit models can perform complex real-time computations on continuous input streams, provided that these computations can be carried out with a rapidly fading memory. We investigate the computational capability of such circuits in the more realistic case where not only readout neurons, but in addition a few neurons within the circuit, have been trained for specific tasks. This is essentially equivalent to the case where the output of trained readout neurons is fed back into the circuit. We show that this new model overcomes the limitation of a rapidly fading memory. In fact, we prove that in the idealized case without noise it can carry out any conceivable digital or analog computation on time-varying inputs. But even with noise, the resulting computational model can perform a large class of biologically relevant real-time computations that require a nonfading memory. We demonstrate these computational implications of feedback both theoretically, and through computer simulations of detailed cortical microcircuit models that are subject to noise and have complex inherent dynamics. We show that the application of simple learning procedures (such as linear regression or perceptron learning) to a few neurons enables such circuits to represent time over behaviorally relevant long time spans, to integrate evidence from incoming spike trains over longer periods of time, and to process new information contained in such spike trains in diverse ways according to the current internal state of the circuit. In particular we show that such generic cortical microcircuits with feedback provide a new model for working memory that is consistent with a large set of biological constraints. Although this article examines primarily the computational role of feedback in circuits of neurons, the mathematical principles on which its analysis is based apply to a variety of dynamical systems. Hence they may also throw new light on the computational role of feedback in other complex biological dynamical systems, such as, for example, genetic regulatory networks.