Insights into the relation between noise and biological complexity
Abstract
Understanding under which conditions the increase of systems complexity is evolutionary advantageous, and how this trend is related to the modulation of the intrinsic noise, are fascinating issues of utmost importance for synthetic and systems biology. To get insights into these matters, we analyzed chemical reaction networks with different topologies and degrees of complexity, interacting or not with the environment. We showed that the global level of fluctuations at the steady state, as measured by the sum of the Fano factors of the number of molecules of all species, is directly related to the topology of the network. For systems with zero deficiency, this sum is constant and equal to the rank of the network. For higher deficiencies, we observed an increase or decrease of the fluctuation levels according to the values of the reaction fluxes that link internal species, multiplied by the associated stoichiometry. We showed that the noise is reduced when the fluxes all flow towards the species of higher complexity, whereas it is amplified when the fluxes are directed towards lower complexity species.
pacs:
02.50.Ey, 05.10.Gg, 05.40.Ca, 87.18.hI Introduction
Fluctuations play a major role in the dynamics of a large variety of complex biological processes. To properly describe the stochastic and heterogeneous nature of systems such as the transcriptional machinery (1) or cell differentiation processes (2), stochastic modeling approaches are indispensable (3). For example, intensive efforts have been devoted in the last decade to the characterization of stochasticity in chemical reaction networks (CRNs) by studying the propagation of fluctuations through the networks (4); (5) or by analyzing the relation between thermodynamic properties and noise levels (6); (7). However, the full understanding of the stochastic properties of biomolecular networks remains an intricate goal that is far from being met.
A challenging open question concerns the relationship between the complexity of biological systems and the modulation of the intrinsic noise. At first glance, the evolutionary pressure, which led from unicellular organisms to higher eukaryotes, tends to favor both complexity increase and noise reduction, but this is clearly not a general rule. On the basis of modeling investigations, it has been shown that, for some systems, the increase of systems complexity is associated to a decrease of the noise level (8), whereas other processes are noisedriven (9); for still other systems, the intrinsic noise can be either increased or decreased according to the parameter values of the model (10).
The lack of a general understanding is due to the fact that CRNs describing biological systems are usually very large and complex and thus complicated to model mathematically, especially via stochastic simulations in which the parameter space becomes rapidly too large to be tractable. Hence, only toy models can be realistically analyzed.
The purpose of this work is to deepen our knowledge on noise modulation in stochastic CRNs. This was done by analyzing model systems with different degrees of complexity and by exploring analytically and numerically their dynamical behavior and parameter spaces, using the Itō stochastic differential equations formalism. For the systems studied, we observed general relations between some structural characteristics of the CRNs and the noise levels evaluated by the Fano factors of the biochemical species involved. We conjecture the validity of these relations for general classes of CRNs.
Ii Stochastic Chemical Reaction Networks
Let us briefly review some basic concepts of chemical reaction networks. They are defined by the triplets consisting of sets of chemical species , complexes and elementary reactions . In the case of open systems, the environment is considered as a complex with vanishing stoichiometric coefficients. Let denote the number of molecules of species (with card()) at time . It satisfies the chemical Langevin equation (CLE) in the Itō formalism (11):
(1) 
where is the stoichiometry matrix, the rate of reaction and statistically independent Gaussian white noises (Wiener processes). This equation describes the temporal evolution of and implies that its conditioned probability density function obeys the associated FokkerPlanck equation. Here we assume a massaction kinetic scheme for the CRNs, which means that the rate of a chemical reaction is proportional to the product of the concentrations of the reactants raised to powers equal to their stoichiometric coefficients.
CRNs are said to be complex balanced (12); (13) if, for each complex , the sum of the mean reaction rates for the reactions for which is a reactant complex is equal to the sum of the mean reaction rates for for which is a product complex at the steady state :
(2) 
Detailed balanced CRNs are a subclass of complex balanced systems for which this relation holds separately for each pair of forward and inverse reactions linking two complexes. Detailed balanced steady states correspond to thermodynamic equilibrium states, whereas the others are nonequilibrium steady states (NESS).
To characterize the topology of a given CRN, the notion of deficiency has been introduced (12); (13) as
(3) 
where is the number of linkage classes, namely the number of connected components of the reaction network, and is the dimension of the stoichiometry subspace, namely the rank of the network.
CRNs are also characterized by their reversibility properties. They are said to be reversible if the existence of a reaction that turns one complex into another implies the existence of the reverse reaction. They are weakly reversible if the existence of a reaction path from complex one to complex two implies the existence of a path from complex two to complex one.
We focused here on open and closed, weakly reversible, massaction CRNs of arbitrary deficiency, which admit a steady state that is not always complex balanced.
Iii Model Systems
We studied three types of weakly reversible CRNs, depicted in Fig.1. The simplest one (Fig.1a) is defined by the reaction chain:
(4) 
It represents for example the assembly of monomeric molecules X into homooligomers Y, where indicates the number of monomers in each oligomer. When , this CRN represents the interconversion between two states of the same molecule (e.g. activated or not) or between two localizations (e.g. intra or extracellular). Both species are linked to the environment.
The system of Itō stochastic differential equations (SDE) that describes this CRN reads as (see Eq. (1)):
(5) 
where and indicate the number of molecules of species X and Y. and are the production terms for the variables and respectively, and are the associated degradation terms, while and are interconversion terms. Assuming massaction kinetics, these terms can be written as:
(6) 
where are the parameters of the model. We assumed here that the number of molecules is large enough so that the approximation holds.
We also studied more complex CRNs, which model for instance a biological system in which a molecular species undergoes a homooligomerization process through an intermediate oligomerization step of lower order (Fig.1b):
(7) 
with , or in which the homooligomerization process occurs with or without an intermediate step (Fig.1c):
(8) 
The SDEs modeling these CRNs can be easily obtained by generalizing Eqs. (56).
Iv Intrinsic Noise
A central parameter that quantifies the role of fluctuations in biochemical systems is the Fano factor , defined as the ratio between the variance of a stochastic variable and its mean: . If the variable follows a Poisson distribution, its Fano factor is equal to one. When is larger than one, the intrinsic noise affects more strongly the variable concentration, and the distribution is called superPoissonian. The distribution is subPoissonian when .
To analyze the role of the fluctuations in the different types of CRNs depicted in Fig.1, we computed the sum of the Fano factors of all species at the steady state as a function of the system’s parameters. This Fano factor sum is taken to represent the global noise level of the system. We obtained its expression from the SDEs of Eqs (56) by employing an EulerMayurama time discretization scheme (14). The equations were studied analytically with a moment closure approximation (15) as well as through numerical simulations. A detailed analysis of these systems will be presented in an upcoming paper (16); here we show the main results.
iv.1 Deficiency zero
We studied a series of CRNs among which the following three classes :
 [label=()]

Closed CRNs described in Figs 1ab with for all species ;

Open CRNs described in Figs 1ab with for all but one species;

Open CRNs described in Figs 1ac with .
Types (a) and (b) are detailed balanced systems, admitting equilibrium states for all reaction rates, while (c) is only complex balanced.
For all these CRNs we found analytically that the sum of the steady state Fano factors over all species is equal to the rank :
(9) 
It can be checked that all CRNs with deficiency zero satisfy this simple relation. It is directly related to the result shown for this kind of CRNs in (17): the steady state probability distribution of the number of molecules of each species can be expressed as a product of Poisson distributions, or as a multinomial distribution in the case some conservation laws hold and thus the state space is reduced, for example in closed systems where the number of molecules is conserved.
iv.2 Deficiency one
For systems with deficiency larger than zero that admit a stationary solution, the situation becomes more complicated due to the correlation between the chemical species in the stationary state.
Consider first the homooligomerization CRN shown in Fig. 1a and Eq. (4) with and the two species X and Y linked to the environment. Note that in this model, both monomers and oligomers can be degraded, but also produced; to obtain the more realistic case in which only the monomer is produced, the oligomer production rate must be set to zero.
Using the same analytical procedure as for the case, we showed that the internal flux that flows between the two species X and Y:
(10) 
is generally nonzero. The other mean fluxes of the system involve exchanges with the environment and are proportional to . We obtained the analytical relation for the sum of the Fano factors in terms of the rank and the mean flux:
(11) 
where is a positive function of the system’s parameters. The rank is here equal to . As expected, this equation reduces to Eq. (9) when (i.e. the system is complex balanced) or (i.e. the steady state is detailed balanced). The analytical and numerical results for the sum of Fano factors as a function of , for different values of and the parameter , are plotted in Fig. 2.
Relation (11) means that there is a global reduction of the intrinsic noise level, as measured by the sum of Fano factors, when the mean flux is positive, thus when the net flux is directed towards the chemical species with higher degree of complexity (here the oligomer Y). When the flux is directed towards the species with lower complexity (here the monomer X), there is a global amplification of the noise. We would like to emphasize that the larger the complexity level, the larger the reduction or amplification effect. Indeed, the proportionality coefficient between the Fano factor sum and the flux increases in absolute value with .
This relation can easily be generalized to the other systems depicted in Fig. 1bc, for the parameter values for which there is only one independent mean flux linking two complexes of different reaction stoichiometry. Note that, in the case of several (dependent) fluxes, the positiveness of the proportionality function is only ensured if they all flow from lower to higher complexity or vice versa.
iv.3 Higher deficiencies
To get insights into the intrinsic noise in more complex CRNs and generalize Eq. (11), we studied networks with higher deficiency values. First we considered the CRN of deficiency shown in Fig. 1b and Eq. (7), with and the three species X, Y and Z linked to the environment. It has two independent fluxes that flow between pairs of species:
(12) 
The fluxes that link the species to the environment are linear combinations of these two fluxes.
We also considered the more complex case with deficiency shown in Fig. 1c and Eq. (8). For and all species linked to the environment, this system has three independent fluxes, the two given by Eq. (12) and a third one:
(13) 
We showed analytically and numerically that the global intrinsic noise in these two CRNs, defined as the sum of the Fano factors over all species, is given by :
(14) 
where all are positive functions of the parameters and where () is the reaction’s stoichiometry, equal to , and for , respectively; the rank is . The sum is over the independent internal mean fluxes expressed as a function of forward and backward reactions as , for which the stoichiometry of reactant and product differ. The analytical results for system (7) are plotted in Fig. 3.
Eq. (14) shows that when all the fluxes flow from the species of lower degree of complexity to the species of higher complexity, the global intrinsic noise level is reduced. In contrast, when the fluxes are directed towards the lower complexity species, it is amplified. Again, the effect on the noise is higher for larger stoichiometry values (). When the fluxes do not all have the same sign, there start to be a competition between the fluxes. In that case, whether the noise is reduced or amplified depends on the parameters of the system, as pictorially shown in Fig. 3.
V Conclusion
In this letter we analyzed the global intrinsic noise in CRNs through the estimation of the sum of Fano factors over all chemical species involved. We found that this quantity depend crucially on the value of the CRN’s deficiency. For all weakly reversible CRNs with = 0 the sum of Fano factors is always constant and equal to the rank of the system independently of the model’s parameter. For higher deficiency systems, additional terms appear which are proportional to the fluxes between the complexes times the reaction stoichiometry. If all fluxes flow in the direction of higher complexity, a global reduction of the noise is observed, while an amplification occurs when fluxes are directed towards lower complexity.
To get insights into the biological meaning of our results, consider for example the system composed of monomeric proteins that undergo an homooligomerization process. In this case, that corresponds to Fig. 1a with , the mean flux is directed towards higher complexity and thus the sum of Fano factors is always smaller than or equal to the rank, which signals global noise reduction. In contrast, for systems modeling the chemical hydrolysis of cellulose and hemicellulose into monomeric sugars, where instead and the flux is directed towards the lowest complexity species, we always have global noise amplification.
Several points remain to be addressed. The analysis of the Fano factors of each individual chemical species in the CRN are left for an upcoming paper (16). We will also extend our results to systems with more general kinetic schemes (18); (19); (20). Last but not least, a clear interpretation of our results in terms of entropy production rates will contribute to deepen our physical understanding of the noise modulation in CRNs (6); (7).
Acknowledgements.
FP is Postdoctoral researcher and MR Research Director at the Belgian Fund for Scientific Research (FNRS).References
 M. Elowitz, A. Levine, E. Siggia, and P. Swain, Stochastic gene expression in a single cell, Science 297, 11831186 (2002).
 G. Balazsi, A. van Oudenaarden, J.J. Collins, Cellular DecisionMaking and Biological Noise: From Microbes to Mammals, Cell, 144, 910925 (2011).
 D.J. Wilkinson, Stochastic modelling for quantitative description of heterogeneous biological systems,Nature Reviews Genetics 10, 122133 (2009).
 DF Anderson, JC Mattingly, HF Nijhout, MC Ree, Propagation of Fluctuations in Biochemical Systems, I: Linear SSC Networks, Bulletin of Mathematical Biology, 69,17911813 (2007).
 DF Anderson, JC Mattingly, Propagation of fluctuations in biochemical systems, II: Nonlinear chains, IET Syst Biol 1,31325 (2007).
 M Polettini, A Wachtel, M Esposito, Dissipation in noisy chemical networks: The role of deficiency, The Journal of Chemical Physics 143, 184103 (2015).
 R Rao, M Esposito, Nonequilibrium Thermodynamics of Chemical Reaction Networks: Wisdom from Stochastic Thermodynamics, Phys Rev X 6, 041064 (2016).
 L Cardelli, A CsikászNagy, N Dalchau, M, M Tschaikowski5, Noise Reduction in Complex Biological Switches, Scientific Reports 6, 20214 (2016).
 M Hoffmann, HH Chang, S Huang, DE Ingber, M, Loeffler, J Galle, Noisedriven stem cell and progenitor population dynamics, PLoSOne 3,e2922 (2008).
 M Rooman, J Albert, M Duerinckx. Stochastic noise reduction upon complexification: Positively correlated birthdeath type systems, Journal of Theoretical Biology 354, 113123 (2014).
 Daniel T. Gillespie,The chemical Langevin equatio, Journal of Chemical Physics, 113, 297306 (2000).
 Martin Feinberg, Complex balancing in general kinetic systems, Arch. Rational Mech. Anal. 49, 187194 (1972).
 Fritz J. M. Horn, Necessary and sufficient conditions for complex balancing in chemical kinetics, Arch. Rat. Mech. Anal. 49, 172186 (1972).
 Kloeden PE, Platen E, Numerical Solution of Stochastic Differential Equations, Springer, Berlin (1992).
 Kuehn C, Moment Closure  A Brief Review, Control of SelfOrganizing Complex Systems, Springer (2016).
 Pucci F, Rooman M, The role of the noise in the evolution of biological complexity, in preparation
 DF Anderson, G Craciun, TG Kurtz, Productform stationary distributions for deficiency zero chemical reaction networks, Bulletin of Mathematical Biology, 72, 19471970 (2010).
 AC Barato, U Seifert, Universal bound on the Fano factor on enzyme kinetic, The Journal of Physical Chemistry B, 119, 65556561 (2015).
 JR Moffitt, C Bustamante, Extracting signal from noise : kinetic mechanisms from a MichaelisMenten like expression for Enzymatic fluctuation, FEBS J, 281, 498 (2014).
 DF Anderson, SL Cotter, ProductForm Stationary Distributions for Deficiency Zero Networks with Nonmass Action Kinetics, Bulletin of Mathematical Biology 78, 23902407 (2016)