The Chemical Reaction Networks (CRN) interpreted through the differential semantics, even when restricted to elementary reactions with mass action law kinetics, form a Turing-complete language. This means that any computable real function can thus be programmed, and in fact compiled, in an abstract CRN that will compute it with an arbitrarily high precision. In this computational framework, the information carriers are the molecular concentrations, the required precision is given as input, and the output concentration is guaranteed to satisfy the required precision. On the other hand, one can be interested in estimating the derivative of an unknown input signal or in reading the concentration value of an input molecular species. By nature, such problems can only be approximated with a finite precision. Hence, the computation framework proposed previously cannot be applied and we need to design and analyze custom CRNs to perform these tasks. In this paper, we present an analog-dyadic converter CRN which takes as input one molecular concentration (in [0, 1] but not necessarily computable), and produces as output a sequence of ”on” and ”off” spikes corresponding to some extent to the sequence of bits in the dyadic representation of the input concentration. We provide a detailed analysis of the source of errors and their behavior when varying the reactions rate constants. We conclude by sketching a possible design for a reader module that takes as input an arbitrary concentration and a desired precision and outputs a dyadic encoding approximating the value of the concentration with the desired precision. We leave as an open question to prove the correctness of our construction.

Chemical Reaction Networks (CRNs) are a standard formalism used in chemistry and biology to model complex molecular interaction systems. In the perspective of systems biology, they are a central tool to analyze the high-level functions of the cell in terms of their low-level molecular interactions. In the perspective of synthetic biology, they constitute a target programming language to implement in chemistry new functions either in vitro, in artificial vesicles, or in living cells. In this paper, we describe the CRN synthesis tool part of our CRN modeling and analysis software BIOCHAM (Biochemical Abstract Machine). This compiler transforms any elementary (resp. algebraic) real function into a formal finite CRN to compute it (resp. with absolute functional robustness), through a pipeline of symbolic computation steps, among which quadratization optimization plays a key role to restrict to elementary reactions with at most two reactants and a minimum number of molecular species.

The Turing completeness of continuous Chemical Reaction Networks (CRNs) states that any computable real function can be computed by a continuous CRN on a finite set of molecular species, possibly restricted to elementary reactions, i.e. with at most two reactants and mass action law kinetics. In this paper, we introduce a more stringent notion of robust online analog computation, and Absolute Functional Robustness (AFR), for the CRNs that stabilize the concentration values of some output species to the result of one function of the input species concentrations, in a perfectly robust manner with respect to perturbations of both intermediate and output species. We prove that the set of real functions stabilized by a CRN with mass action law kinetics is precisely the set of real algebraic functions. Based on this result, we present a compiler which takes as input any algebraic function (defined by one polynomial and one point for selecting one branch of the algebraic curve defined by the polynomial) and generates an abstract CRN to stabilize it. Furthermore, we provide error bounds to estimate and control the error of an unperturbed system, under the assumption that the environment inputs are driven by k-Lipschitz functions.

The online estimation of the derivative of an input signal is widespread in control theory and engineering. In the realm of chemical reaction networks (CRN), this raises however a number of specific issues on the different ways to achieve it. A CRN pattern for implementing a derivative block has already been proposed for the PID control of biochemical processes, and proved correct using Tikhonov’s limit theorem. In this paper, we give a detailed mathematical analysis of that CRN, thus clarifying the computed quantity and quantifying the error done as a function of the reaction kinetic parameters. In a synthetic biology perspective, we show how this can be used to compute online functions with CRNs augmented with an error correcting delay for derivatives. In the systems biology perspective, we give the list of models in BioModels containing (in the sense of subgraph epimorphisms) the core derivative CRN, most of which being models of oscillators and control systems in the cell, and discuss in detail two such examples: one model of the circadian clock and one model of a bistable switch.

The Turing completeness of continuous chemical reaction networks (CRNs) states that any computable real function can be computed by a continuous CRN on a finite set of molecular species, possibly restricted to elementary reactions, i.e. with at most two reactants and mass action law kinetics. In this paper, we introduce a notion of online analog computation for the CRNs that stabilize the concentration of their output species to the result of some function of the concentration values of their input species, whatever changes are operated on the inputs during the computation. We prove that the set of real functions stabilized by a CRN with mass action law kinetics is precisely the set of real algebraic functions.

In this short paper extracted from [7], we present a polynomialization algorithm of quadratic time complexity to transform a system of elementary differential equations in polynomial differential equations (PODE). This algorithm is used as a front-end transformation in a pipeline to compile any elementary mathematical function, either of time or of some input variable, into a finite Chemical Reaction Network (CRN) which computes it. We illustrate the performance of our compiler on a benchmark of elementary functions which serve as formal specification of CRN design problems in synthetic biology, and as comparison basis with natural CRNs exhibiting similar behaviours.

The Turing completeness result for continuous chemical reaction networks (CRN) shows that any computable function over the real numbers can be computed by a CRN over a finite set of formal molecular species using at most bimolecular reactions with mass action law kinetics. The proof uses a previous result of Turing completeness for functions defined by polynomial ordinary differential equations (PODE), the dual-rail encoding of real variables by the difference of concentration between two molecular species, and a back-end quadratization transformation to restrict to elementary reactions with at most two reactants. In this paper, we present a polynomialization algorithm of quadratic time complexity to transform a system of elementary differential equations in PODE. This algorithm is used as a front-end transformation to compile any elementary mathematical function, either of time or of some input species, into a finite CRN. We illustrate the performance of our compiler on a benchmark of elementary functions relevant to CRN design problems in synthetic biology specified by mathematical functions. In particular, the abstract CRN obtained by compilation of the Hill function of order 5 is compared to the natural CRN structure of MAPK signalling networks.

Chemical reaction networks (CRNs) are a standard formalism used in chemistry and biology to reason about the dynamics of molecular interaction networks. In their interpretation by ordinary differential equations, CRNs provide a Turing-complete model of analog computattion, in the sense that any computable function over the reals can be computed by a finite number of molecular species with a continuous CRN which approximates the result of that function in one of its components in arbitrary precision. The proof of that result is based on a previous result of Bournez et al. on the Turing-completeness of polyno-mial ordinary differential equations with polynomial initial conditions (PIVP). It uses an encoding of real variables by two non-negative variables for concentrations, and a transformation to an equivalent quadratic PIVP (i.e. with degrees at most 2) for restricting ourselves to at most bimolecular reactions. In this paper, we study the theoretical and practical complexities of the quadratic transformation. We show that both problems of minimizing either the number of variables (i.e., molecular species) or the number of monomials (i.e. elementary reactions) in a quadratic transformation of a PIVP are NP-hard. We present an encoding of those problems in MAX-SAT and show the practical complexity of this algorithm on a benchmark of quadratization problems inspired from CRN design problems.

Numerous biological systems are known to harbor a form of logarithmic behavior, from Weber's law to bacterial chemotaxis. Such a log-response allows for sensitivity to small relative variations of biochemical inputs over a large range of concentration values. Here we use a genetic algorithm to evolve biochemical networks displaying a logarithmic response. A quasi-perfect log-response implemented by the same core network evolves in a convergent way across our different in silico replications. The best network is able to fit a logarithm over 4 orders of magnitude with an accuracy of the order of 1%. At the heart of this network, we show that a logarithmic approximation may be implemented with one single nonlinear interaction, that can be interpreted either as multisite phosphorylations or as a ligand induced multimerization. We provide an analytical explanation for the effect and exhibit constraints on parameters. Biological log-response might thus be easier to implement than usually assumed.

One goal of synthetic biology is to implement useful functions with biochemical reactions, either by reprogramming living cells or programming artificial vesicles. In this perspective, we consider Chemical Reaction Networks (CRN) as a programming language, and investigate the CRN program synthesis problem. Recent work has shown that CRN interpreted by differential equations are Turing-complete and can be seen as analog computers where the molecular concentrations play the role of information carriers. Any real function that is computable by a Turing machine in arbitrary precision can thus be computed by a CRN over a finite set of molecular species. The proof of this result gives a numerical method to generate a finite CRN for implementing a real function presented as the solution of a Polynomial Initial Values Problem (PIVP). In this paper, we study an alternative method based on artificial evolution to build a CRN that approximates a real function given on finite sets of input values. We present a nested search algorithm that evolves the structure of the CRN and optimizes the kinetic parameters at each generation. We evaluate this algorithm on the Heaviside and Cosine functions both as functions of time and functions of input molecular species. We then compare the CRN obtained by artificial evolution both to the CRN generated by the numerical method from a PIVP definition of the function, and to the natural CRN found in the BioModels repository for switches and oscillators.

Molecular networks are at the core of most cellular decisions, but are often difficult to comprehend. Reverse engineering of network architecture from their functions has proved fruitful to classify and predict the structure and function of molecular networks, suggesting new experimental tests and biological predictions. We present φ-evo, an open-source program to evolve in silico phenotypic networks performing a given biological function. We include implementations for evolution of biochemical adaptation, adaptive sorting for immune recognition, metazoan development (somitogenesis, hox patterning), as well as Pareto evolution. We detail the program architecture based on C, Python 3, and a Jupyter interface for project configuration and network analysis. We illustrate the predictive power of φ-evo by first recovering the asymmetrical structure of the lac operon regulation from an objective function with symmetrical constraints. Second, we use the problem of hox-like embryonic patterning to show how a single effective fitness can emerge from multi-objective (Pareto) evolution. φ-evo provides an efficient approach and user-friendly interface for the phenotypic prediction of networks and the numerical study of evolution itself.

We consider the general problem of sensitive and specific discrimination between biochemical species. An important instance is immune discrimination between self and not-self, where it is also observed experimentally that ligands just below the discrimination threshold negatively impact response, a phenomenon called antagonism. We characterize mathematically the generic properties of such discrimination, first relating it to biochemical adaptation. Then, based on basic biochemical rules, we establish that, surprisingly, antagonism is a generic consequence of any strictly specific discrimination made independently from ligand concentration. Thus antagonism constitutes a 'phenotypic spandrel': a phenotype existing as a necessary by-product of another phenotype. We exhibit a simple analytic model of discrimination displaying antagonism, where antagonism strength is linear in distance from the detection threshold. This contrasts with traditional proofreading based models where antagonism vanishes far from threshold and thus displays an inverted hierarchy of antagonism compared to simpler models. The phenotypic spandrel studied here is expected to structure many decision pathways such as immune detection mediated by TCRs and FCϵRIs, as well as endocrine signalling/disruption.

The sequence of a protein is not only constrained by its physical and biochemical properties under current selection, but also by features of its past evolutionary history. Understanding the extent and the form that these evolutionary constraints may take is important to interpret the information in protein sequences. To study this problem, we introduce a simple but physical model of protein evolution where selection targets allostery, the functional coupling of distal sites on protein surfaces. This model shows how the geometrical organization of couplings between amino acids within a protein structure can depend crucially on its evolutionary history. In particular, two scenarios are found to generate a spatial concentration of functional constraints: high mutation rates and fluctuating selective pressures. This second scenario offers a plausible explanation for the high tolerance of natural proteins to mutations and for the spatial organization of their least tolerant amino acids, as revealed by sequence analysis and mutagenesis experiments. It also implies a faculty to adapt to new selective pressures that is consistent with observations. The model illustrates how several independent functional modules may emerge within the same protein structure, depending on the nature of past environmental fluctuations. Our model thus relates the evolutionary history of proteins to the geometry of their functional constraints, with implications for decoding and engineering protein sequences.