Exponential equilibration of genetic circuits using entropy methods

Journal of Mathematical Biology (2019) 78:373–411 https://doi.org/10.1007/s00285-018-1277-z Mathematical Biology Exponential equilibration of genetic circuits using entropy methods José A. Cañizo1·José A. Carrillo2·Manuel Pájaro3 Received: 24 January 2018 / Revised: 16 July 2018 / Published online: 17 August 2018 © The Author(s) 2018 Abstract We analyse a continuum model for genetic circuits based on a partial integrodifferential equation initially proposed in Friedman et al. (Phys Rev Lett 97(16):168302, 2006) as an approximation of a chemical master equation. We use entropy methods to show exponentially fast convergence to equilibrium for this model with explicit bounds. The asymptotic equilibration for the multidimensional case of more than one gene is also obtained under suitable assumptions on the equilibrium stationary states. The asymptotic equilibration property for networks involving one and more than one gene is investigated via numerical simulations. Mathematics Subject Classification 35B40 ·92Dxx ·39B99 ·65M99 1 Introduction Translation of the information encoded in genes is responsible for all cellular functions. The decoding of DNA can be summarised, following the central dogma of molecular biology, in two steps: the transcription into messenger RNA and the translation into proteins. Cells produce responses to environmental signals, thanks to the regulation of DNA expression via certain feedback mechanism activating or inhibiting the genes. Typically, regulation is produced by the union of proteins to the DNA binding sites. BJosé A. Carrillo [email protected] José A. Cañizo [email protected] Manuel Pájaro [email protected] 1Departamento de Matemática Aplicada, Universidad de Granada, 18071 Granada, Spain 2Department of Mathematics, Imperial College London, London, SW7 2AZ, UK 3BioProcess Engineering Group, IIM-CSIC, Spanish Council for Scientific Research, Eduardo Cabello 6, 36208 Vigo, Spain 123 374 J. A. Cañizo et al. Moreover, the number of species involved in gene regulatory networks (gene expression together with their regulation) is small, which makes its behaviour inherently stochastic (Elowitz et al. 2002; Gillespie 2007; Kaeet al. 2005; McAdams and Arkin 1997; Paulsson 2004). This underlying stochastic behaviour in gene regulatory networks is captured by using the chemical master equation (CME) (Kepler and Elston 2001; Mackey et al. 2011; Paulsson 2005; Sherman and Cohen 2014). However, the CME solution is unavailable in most cases, due to the large (even infinite) number of coupled equations. There are two main ways to obtain the CME solution: via stochastic simulation or via approximations of the CME. One of the most extended methods to reproduce the CME dynamics using stochastic realisations is the stochastic simulation algorithm (SSA) (Gillespie 1976,2007). This method has no restrictions in its applicability, even though it is computationally expensive. On the other hand, CME approximations which remain valid under certain conditions include the finite state projection (Munsky and Khammash 2006), moment methods (Engblom 2006; Hasenauer et al. 2015), linear noise approximations (Thomas et al. 2014; Kampen 2007; Wallace et al. 2012)or hybrid models (Jahnke 2011). In addition to the above mentioned methods, assuming that protein production takes place in bursts one can obtain a partial integro-differential equation (PIDE) as a continuous approximation of the CME. This PIDE has a mathematical structure very similar to kinetic and transport equations in mathematical biology (Perthame 2007) and it admits an analytical solution for its steady state in the case of networks involving only one gene. In the next subsections, we describe both the one dimensional PIDE model (Friedman et al. 2006) for self-regulated gene networks and the generalised PIDE model (Pájaro et al. 2017) for arbitrary genetic circuits. We will discuss the main properties of the stationary states in one dimension to finally explain the main results of this work. 1.1 1-dimensional PIDE model The kinetic equation, first proposed by Friedman et al. (2006), is a continuous approximation of the CME for gene self-regulatory networks. A schematic representation of this genetic circuit is illustrated in Fig. 1, where the transcription-translation mechanism from DNA to a protein Xis shown. Note that DNA transcribes into messenger RNA not only from the active state at rate (per unit time τ)km, but also from the inactive state with rate constant kεlower than km, which is known as basal transcription level or transcriptional leakage (Friedman et al. 2006; Ochab-Marcinek and Tabaka 2015; Pájaro et al. 2015). The messenger RNA transcribes into protein Xfollowing a first-order process with rate constant (per unit time) kx. The messenger RNA and protein are degraded at rate constants γmand γxrespectively. For self-regulated gene networks, activation or inhibition of the DNA promoter is produced by the union of the protein expressed to the DNA binding sites (feedback mechanism). So that, under protein action the promoter can switch between its inactive (DNAoff ) and active (DNAon) forms, with rate constants kon and koff respectively (see Fig. 1). There are two types of feedback mechanism: positive or negative, cor123 Exponential equilibration of genetic circuits using entropy methods 375 DNAoff k kon koff DNAon kmmRNA γm kxX γx X∅∅ Fig. 1 Schematic representation of the transcription-translation mechanism under study. The promoters associated with the gene of interest are assumed to switch between active (DNAon) and inactive (DNAoff ) states, with rate constants kon and koff per unit time, respectively. In this study, the transition is assumed to be controlled by a feedback mechanism induced by the binding/unbinding of a given number of X-protein molecules, what makes the network self-regulated. Transcription of messenger RNA (mRNA) from the active DNA form, and translation into protein Xare assumed to occur at rates (per unit time) kmand kx, respectively. kεis the rate constant associated with transcriptional leakage. The mRNA and protein degradations are assumed to occur by first order processes with rate constants γmand γx, respectively responding to whether the protein inhibits or promotes their production, respectively. The fraction of the promoter in the active or inactive state is typically described by Hill functions (Alon 2007). We can express the probability that the promoter is in its inactive state as a function of the protein amount x, denoted by ρ:R+→[0,1](see Ochab-Marcinek and Tabaka 2015; Pájaro et al. 2015): ρ(x)=xH xH+KH,(1.1) where K:= koff kon is the equilibrium binding constant and H∈Z\{0}is the Hill coefficient which is positive if Hproteins bound to the DNA inhibiting their production (negative feedback) and negative if |H|proteins bound to the DNA activating their production (positive feedback). Then, the rate RTof messenger RNA production (transcription) can be written as function of the Hill expression (1.1), RT=kmc(x), with the input function c(x):= (1−ρ(x))+ρ(x)ε, where εis the leakage constant defined as ε:= kε km. Note that the function RTaccounts for the messenger RNA production both from the DNA active state (with probability 1 −ρ(x)) with rate constant kmand from the inactive DNA (with probability ρ(x)) with lower rate constant kε. The PIDE model is valid under the assumption of protein production in bursts. So, we consider gene self-regulatory networks where the degradation rate of mRNA is much faster than the corresponding to protein, γm/γx1. Such condition is verified in many gene regulatory networks, both in prokaryotic and eukaryotic organisms (Shahrezaei and Swain 2008; Dar et al. 2012), and results in protein being produced in bursts. As suggested in Friedman et al. (2006) and Elgart et al. (2011), the burst size (denoted by b=kx γm) is typically modelled by an exponential distribution. The conditional probability for protein level to jump from a state yto a state x>yafter a burst is proportional to: ω(x−y)=1 bexp −x−y b,for x>y>0.(1.2) 123 376 J. A. Cañizo et al. The temporal evolution of the probability density function of the amount of proteins, p:R+×R+→R+is described by the following PIDE model: ∂p ∂t(t,x)−∂(xp) ∂x(t,x)=ax 0 ω(x−y)c(y)p(t,y)dy−ac(x)p(t,x), (1.3) where τis time, t=γxτrepresents a dimensionless time associated to the time scale of protein degradation, a=km γxis the dimensionless rate constant related to transcription, which represents the mean number of bursts (burst frequency) and ω(x−y)is given by (1.2). The input function c:R+→[ε, 1], which represents the feedback mechanism, takes the form (Ochab-Marcinek and Tabaka 2015; Pájaro et al. 2015): c(x)=KH+εxH KH+xH,x>0.(1.4) Note that the above input function can be constant, equal to one, when the protein does not promote or repress its production (open loop). This constant c(x)=1isused when the DNA is always in its active state, thus implying a unique messenger RNA production rate (km), reducing the system complexity. We denote the stationary solution of Eq. (1.3) (which we sometimes call equilibrium)asP∞(x), which therefore verifies the following equation: ∂[xP ∞(x)] ∂x=−ax 0 ω(x−y)c(y)P∞(y)dy+ac(x)P∞(x). (1.5) We say a stationary solution is normalised when its integral over [0,+∞)(which we sometimes call its mass) is equal to 1. This equation has a unique solution with mass 1, which can be written out explicitly as (Ochab-Marcinek and Tabaka 2015; Pájaro et al. 2015): P∞(x):= Z[ρ(x)]a(1−ε) Hx−(1−aε)e−x b=ZxH+KHa(ε−1) Hxa−1e−x b,(1.6) with ρ(x)defined in (1.1) and Zbeing a normalising constant such that ∞ 0P∞(x)dx= 1. Alternatively, stationary solutions may be studied by considering the zero-flux case; see for example Bokes and Singh (May 2017); Bokes et al. (Jul 2018). In case of no self-regulation (open loop network with c(x)=1; that is, =1) the stationary solution is a gamma distribution (Friedman et al. 2006), which is in fact the limit of (1.6)astends to 1: P∞(x):= xa−1e−x/b ba(a).(1.7) 123 Exponential equilibration of genetic circuits using entropy methods 377 1.2 Generalised n-dimensional PIDE model Recently the 1D PIDE model has been extended to overcome more general gene regulatory networks than the self-regulation considered by Friedman et al. (2006). As a first step in this extension, Bokes and Singh (2015) propose the use of variable protein degradation rate, in order to accommodate gene networks with decoy binding sites (Lee and Maheshri 2012) to the PIDE model structure. Finally, including the previous models and considering genetic networks involving more than one gene Pájaro et al. (2017) proposed the generalised PIDE model for any number of genes. In Pájaro et al. (2017) a general gene regulatory network comprising ngenes, G={DNA1,...,DNAi,...,DNAn}, is proposed. These genes encoded by DNA-subchains are transcribed into ndifferent messenger RNAs M={mRNA1, ...,mRNAi,...,mRNAn}, which are translated into nproteins types X={X1,..., Xi,...,Xn}. We show a schematic representation of the general network in Fig. 2, which is similar to the self-regulation circuit. The main differences are that: (i) each DNA type can be regulated by others different proteins than the one expressed by the considered gene (cross regulation), and (ii) the protein degradation rate can be a variable function of all proteins types considered. The structure of this multidimensional network is equivalent to the previous selfregulation case. Each promoter can switch from the inactive states (DNAioff )tothe active one (DNAion) or vice versa with rate constants ki on and ki off respectively. The leakage (basal) messenger RNA production from the inactive promoter is conserved at lower rate constant (ki ε) than its production from the active state (ki m). Each imessenger RNA type is translated into the protein Xiat rate constant ki x. Both messengers RNA and proteins are degraded with rates γi mand γi x(x)respectively. Note that for this general network the total rate of production of mRNAi,Ri T, can be written as the rate constant production from the active DNAistate times one input function ci(x)describing all possible types of feedback mechanism. However, there are not universal expressions for ci(x), due to their dependence on the regulatory mechanism considered (the messenger RNA production can occur from intermediate DNA states between the total activated and the total repressed ones), some examples have been described in Alon (2007) and Pájaro et al. (2017). Without loss of generality, we can construct the input function verifying that its image is a positive interval, ci:Rn +→[εi,1], where the leakage constant εiis defined as ki ε/ki mwith ki εbeing the mRNAi rate constant from the total repressed DNAi(the lowest rate of mRNAi production). Considering the set of nproteins X={X1,...,Xn}, we define the n-vector x=(x1,...,xn)∈Rn +as the amount of each protein type. The generalised (ndimensional) PIDE model, proposed in Pájaro et al. (2017), describes the temporal evolution of the joint density distribution function of nproteins p:R+×Rn +→R+: ∂p ∂t(t,x)= n  i=1∂ ∂xiγi x(x)xip(x) +ki mxi 0 ωi(xi−yi)ci(yi)p(t,yi)dyi−ki mci(x)p(x)(1.8) 123 378 J. A. Cañizo et al. DNAioff ki ki on ki off DNAion ki mmRN Ai γi m ki xXi γi x(x) XJ∅∅ Fig. 2 Schematic representation of the transcription-translation mechanism under study. The promoters associated with the genes of interest are assumed to switch between active (DNAion) and inactive (DNAioff ) states, with rate constants ki on and ki off per unit time, respectively. The transition is assumed to be controlled by a feedback mechanism induced by the binding/unbinding of a given number of Xjprotein molecules with j∈J(more than one protein type can bind to the DNA), which makes the network self-regulated if i=jor cross-regulated if j= i. Transcription of messenger RNA (mRNAi) from the active DNAi form, and translation into protein Xiare assumed to occur at rates (per unit time) ki mand ki x, respectively. ki εis the rate constant associated with transcriptional leakage. The mRNAidegradation is assumed to occur by first order processes with rate constant γi m. Degradation of the Xi-protein may follow different pathways, which is modelled by the function γi x(x), with γi x:Rn +→R+ where yirepresents the vector state xwith its i-th position changed to yi, (that is: (yi)j=xjif j= iand (yi)j=yiif j=i), and γi x(x)is the degradation rate function of each protein. The first term in the right-hand side of the equation accounts for protein degradation whereas the integral describes protein production by bursts. The burst size is assumed to follow an exponential distribution, what leads to the conditional probability for protein jumping from a state yito a state xiafter a burst be given by: ωi(xi−yi)=1 bi exp −xi−yi bi where bi=ki x γi m are dimensionless frequencies associated to translation which corresponds with the mean protein produced per burst (burst size). The function ci(x) (ci:Rn +→[εi,1]) is an input function, which models the regulation mechanism of the network considered. The stationary solution P∞(x)of (1.8) satisfies: n  i=1∂ ∂xiγi x(x)xiP∞(x) +ki mxi 0 ωi(xi−yi)ci(yi)P∞(yi)dyi−ki mci(x)P∞(x)=0.(1.9) Note that an analytical expression for the steady state solution is not known for the general case of the PIDE model (1.8). Some properties of the 1D solution remain valid for the nD steady state since P∞(x)is a probability density function, then Rn +P∞(x)dx=1. However, we do not have any other prior information about the properties of stationary solutions. 123 Exponential equilibration of genetic circuits using entropy methods 379 1.3 Main results In this work we will apply entropy methods in order to analyse the asymptotic equilibration for the kinetic equations (1.3) and (1.8). These equations bear asimilar structure to the self-similar fragmentation and the growth-fragmentation equations (Perthame and Ryzhik 2005; Laurençot and Perthame 2009; Doumic 2010; Cáceres et al. 2011; Balagué et al. 2013), used for instance in cell division modelling. In those cases, the transport term makes the cluster size of particles grow while the integral term breaks the particles into pieces of smaller size. In our present models, the transport term degrades the number density of proteins while the integral term makes the protein number density to grow. In fact, the kinetic equations (1.3) and (1.8) have the structure of linear population models as in Michel et al. (2004,2005) and Carrillo et al. (2011) for which the socalled general relative entropy applies. This fact already reported in Pájaro et al. (2016) implies the existence of infinitely many Lyapunov functionals for these models useful for different purposes among which to analyse their asymptotic behavior. We will make a summary of the main properties of Eq. (1.3) in Sect. 2together with a quick treatment of the well-posedness theory for these models. They are easily generalisable to the multidimensional case (1.8). In Sects. 3and 4, we will improve over the direct application of the general relative entropy method in Pájaro et al. (2016). On one hand, we study in Sect. 3the case of gene circuits involving one gene, Eq. (1.3), a direct functional inequality between the L2-relative entropy and its production leading to exponential convergence. In order to fix our setting, we recall that ωis given by (1.2)forsomeb>0, and c=c(x)is given by (1.4), for some constants K>0, H∈Z\{0}and 0 <≤1; and a>0 is a constant. For 1 ≤p<+∞ we denote by Lp() the usual Lebesgue spaces of real functions fon such that |f|pis integrable in the Lebesgue sense. We also write Lp(, w) to denote the corresponding spaces of functions fsuch that |f|pis integrable with a weight w. Theorem 1.1 (Long-time behaviour for the 1-dimensional model) Let p0be a probability distribution such that p0∈L1((0,+∞)) ∩L2((0,+∞), P−1 ∞), and let p be the mild solution to Eq. (1.3)with initial data p0(see Definition 2.1). There exists a constant λ>0depending only on the parameters of the equation (and not on p0) such that p(t,·)−P∞L2((0,+∞),P−1 ∞)≤e−λtp0−P∞L2((0,+∞),P−1 ∞). The value of λcan be estimated explicitly from the arguments in the proof, though we do not consider the specific value to be a good approximation of the optimal decay rate. The behaviour of the stationary solutions P∞(x)near the origin and infinity is crucial for direct functional inequalities involving the relative entropy and its production in the one dimensional case. What we are showing is essentially a spectral gap in a weighted L2norm, and some remarks are in order regarding the specific choice of space L2((0,+∞), P−1 ∞)that we have made. As will be seen later, this space is very natural for the technique we 123 380 J. A. Cañizo et al. are going to use, since the evolution operator is contractive in this norm, and a similar observation is true for any Markov semigroup with an equilibrium. However, it is very likely that this operator also has a spectral gap in L2norms with different weights, in weighted L1norms, and in other metrics, as is often the case with Markov operators. In many examples (such as the Fokker-Planck equation) it is known that the spectral gap property breaks for weights which are slowly decaying, so that there may not be a spectral gap in L1, for example. In those cases there are well-known examples of initial data with slowly-decaying tails whose associated solution converges to equilibrium as slowly as one wishes. The same happens for example to the Boltzmann equation from kinetic theory; we refer to Gualdani et al. (2010) for details on the extension of spectral gaps to different weights. So the weight is not only a technical assumption: there may be norms and weights in which the convergence is not exponential. However, exponential weights as the ones we use are probably far from being the optimal ones where one can show a similar result. Section 4is devoted to the analysis of the multidimensional Eq. (1.8) corresponding to multiple genes involved in the gene transcription. In this case, solutions to the stationary problem (1.9) are not explicit and hence we are not able to control precisely the behaviour of the stationary solutions near the origin and infinity as before. For this reason, we are only able to show convergence towards a unique equilibrium solution assuming its existence with suitable behavior near the origin and infinity: Theorem 1.2 (Long-time behaviour for the nDmodel) Given any mild solution p with normalised nonnegative initial data p0∈L1(R+)to Eq. (1.8)and given a normalised stationary solution P∞(x)to (1.8)satisfying the technical Assumption 4.1 from Sect. 4, it holds that lim t→∞ Rn + |p(t,x)−P∞(x)|2dx=0. As a consequence, if a normalised stationary solution P∞(x)of (1.8)and satisfying Assumption 4.1 exists, it is unique. The proof is based on a weaker variant of our one-dimensional inequality, in which the control between the relative entropy and its production is obtained except for an error term which happens to be small under the assumptions of the behavior of the stationary solution P∞(x). Both results of equilibration are illustrated with numerical simulations in their corresponding sections. 2 Mathematical preliminaries and entropy methods 2.1 Properties of stationary solutions Let us start by discussing the basic properties of the one dimensional stationary states to (1.3). The behaviour of the stationary state at zero and at +∞ depends on both r=aε−1 and adue to the presence of the function ρ(x)and its dependence on H. It is as follows: 123 Exponential equilibration of genetic circuits using entropy methods 381 b a Binary region Bimodal region 1/ε 12 34 5 Fig. 3 Regions in the parameter space, where protein distribution exhibits different behaviours for H<0. There are two large areas where the protein distribution change fundamentally its properties, the first including the shapes one and two, where a<1 εand limx→0P∞(x)=+∞and the second with P∞(x) finite for all non-negative x, which includes shapes three to five 1. If H>0, then P∞(x)≃xa−1as x→0+and P∞(x)≃xre−x/bas x→+∞. Then the stationary state P∞(x)exhibits a singularity at zero for 0 <a<1 and it is smooth otherwise having zero limit for a>1 and a positive limit for a=1. 2. If H<0, then P∞(x)≃xras x→0+and P∞(x)≃xa−1e−x/bas x→+∞. Then the stationary state P∞(x)exhibits a singularity at zero for aε<1 and it is smooth otherwise having zero limit for aε>1 and a positive limit for aε=1. As a particular case, if c(x)≡1 then P∞(x)is given by (1.7) and we have P∞(x)≃ xa−1as x→0+and P∞(x)≃xa−1e−x/bas x→+∞. Then the stationary state P∞(x)exhibits a singularity at zero for a<1 and it is smooth otherwise having zero limit for a>1 and a positive limit for a=1. Note that in all cases limx→∞ P∞(x)=0. As we can see in Fig 3, the stationary solution has five different qualitative behaviours for H<0(seealsoPájaroetal. 2015): 1. If a<1 ε, then limx→0P∞(x)=∞. 1.1 Only one peak in x=0(Case1Fig 3). 1.2 Two peaks one in x=0 and another in x>0(Case2Fig 3). 2. If a>1 ε, then limx→0P∞(x)=0. If a≥1 ε, then limx→0P∞(x)=Mwith M≥0. 2.1 Only one peak in x>0 but close to x=0(Case3Fig 3). 2.2 Two different peaks at two points x1,x2>0(Case4Fig3). 2.3 Only one peak in x≥0(Case5Fig 3). 123 388 J. A. Cañizo et al. 3 Exponential convergence for the 1D PIDE model In this section our aim is to prove that Eq. (1.3) converges exponentially to the steady state, P∞. For this purpose, we consider the L2-relative entropy, i.e., the convex function His chosen as H(u)=(u−1)2, and G2(u)(t):=∞ 0 P∞(x)(u(t,x)−1)2dx =∞ 0 p2(t,x) P2 ∞(x)P∞(x)dx−1=∞ 0 u2(t,x)P∞(x)dx−1, where we have used that p(t,x)and P∞(x)are probability density functions. Now, by replacing the value of the considered convex function in Proposition 2.4, we obtain the following identity D2(u)(t):= −dG2(u) dt =a∞ 0∞ y ω(x−y)(u(t,x)−u(t,y))2c(y)P∞(y)dxdy.(3.1) The entropy method consists in finding conditions under which the following functional inequality holds: G2(u)≤1 2βD2(u). (3.2) Notice that the dependence on the time variable can be forgotten at this point, since our objective is to show such an inequality among a subset of suitable probability densities. For this purpose, we start by rewriting G2(u)in a equivalent form (Cáceres et al. 2011): Lemma 3.1 Given a non-negative measurable function P∞:(0,∞)→R+such that ∞ 0P∞(x)dx=1and defining the functional H2(u):= ∞ 0∞ y P∞(x)P∞(y)(u(x)−u(y))2dxdy, there holds G2(u)=H2(u). Proof Expanding the square implies G2(u)=∞ 0 P∞(x)(u(x)−1)2dx=∞ 0 P∞(x)u(x)2dx−1,(3.3) while H2(u)is a symmetric function, so that: H2(u)(τ) =1 2∞ 0∞ 0 P∞(x)P∞(y)(u(x)−u(y))2dxdy =1 2∞ 0∞ 0 P∞(x)P∞(y)u(x)2−2u(x)u(y)+u(y)2dxdy 123 Exponential equilibration of genetic circuits using entropy methods 389 =∞ 0∞ 0 P∞(x)P∞(y)u(x)2dxdy −∞ 0∞ 0 P∞(x)P∞(y)u(x)u(y)dxdy =∞ 0 P∞(x)u(x)2∞ 0 P∞(y)dydx−∞ 0∞ 0 p(x)p(y)dxdy =∞ 0 P∞(x)u(x)2dx−1, which is equal to (3.3).  As consequence of this lemma we are reduced to show the inequality H2(u)≤1 2βD2(u), (3.4) among a suitable subset of probability densities. 3.1 Entropy-entropy production inequality We start by obtaining bounds for the steady state solution P∞, of the Friedman Eq. (1.3). Lemma 3.2 (P∞bounds) For δ>0we define the intervals of length 1 2: Ik,δ := δ+k 2,δ+k+1 2,k≥0integer, and pk:= Cδ+k 2H +KHa(ε−1) Hδ+k 2a−1 e −(δ+k 2) b=P∞δ+k 2. Then, the following inequality holds: A(δ) ≤P∞(x) pk ≤B(δ), ∀x∈Ik,δ and ∀k,(3.5) with P∞(x)given by (1.6)and A(δ) and B(δ) being positive constants that only depend on δ(and network parameters), but they are independent of protein amount k. Proof Note that xH+KHa(ε−1) Hand e−x bare decreasing functions, so that their maxima are at ¯x0=δ+k 2and their minima are at ¯x1=δ+k+1 2in Ik,δ.Thetermxa−1 shows different behaviours which depend on the parameter a, (this term is increasing 123 390 J. A. Cañizo et al. if a>1, constant if a=1 and decreasing if a<1). So that, we can bound P∞(x)in the interval Ik,δ as follows: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ g(¯x1)(δ +k 2)a−1≤P∞(x)≤g(¯x0)δ+k+1 2a−1if a>1 g(¯x1)≤P∞(x)≤g(¯x0)if a=1 g(¯x1)δ+k+1 2a−1≤P∞(x)≤g(¯x0)δ+k 2a−1if a<1 (3.6) where g(x)=ZxH+KHa(ε−1) He−x b. Now, in order to calculate the bounds of P∞(x) pk, we divide the expression (3.6)by pkto obtain A(δ, k)≤P∞(x) pk≤B(δ, k)with the functions Aand Bbeing, A(δ, k):= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (δ +k+1 2)H+KH (δ +k 2)H+KHa(ε−1) H e−1 2bif a≥1 (δ +k+1 2)H+KH (δ +k 2)H+KHa(ε−1) H e−1 2b2δ+k+1 2δ+ka−1 if a<1 and B(δ, k):= ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 2δ+k+1 2δ+ka−1 if a>1 1ifa≤1 Notice that, lim k→∞ A(δ, k)=e−1 2b,lim k→∞ B(δ, k)=1, implies that A(δ) := min k≥0(A(δ, k))and B(δ) := max k≥0(B(δ, k))are well-defined and positive, leading to desired inequality (3.5).  Note that inequality (3.5) can be directly checked for the simplest open loop case, whose stationary solution is given by (1.7). Lemma 3.3 Let us define Mj:= j−1  k=1 1 mk ,(3.7) 123 Exponential equilibration of genetic circuits using entropy methods 391 with {mk}k≥1a positive sequence given by mk=pke δ+k 2 2b. Then, there exists C >0 such that mk ∞  j=k+1 Mjpj≤Cpk,for all k ∈N.(3.8) Proof We define {aj}j≥1with aj=1 mjto calculate the following limit lim j→∞ aj+1−aj Mj+1−Mj =lim j→∞ ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ δ+j+1 2H +KHa(1−ε) Hδ+j+1 21−a δ+j 2H +KHa(1−ε) Hδ+j 21−a e1 4b−1⎞ ⎟ ⎟ ⎟ ⎟ ⎠ =e1 4b−1. Since this limit exists and {Mj}j≥1is a strictly increasing and divergent sequence, we can use the Stolz-Cesàro theorem to obtain that Mj≤C0aj, with C0>0 constant. Then, mk ∞  j=k+1 Mjpj≤C0mk ∞  j=k+1 ajpj. The summation term at the right hand side can be calculated as follows ∞  j=k+1 ajpj= ∞  j=k+1 e−2δ+j 4b=e−2b−1 4b e−1e−2δ+k 4b, so that mk ∞  j=k+1 Mjpj≤Cmke−2δ+k 4b=Cpk, with C=C0e−2b−1 4b e−1, concluding the proof.  In order to prove the exponential convergence of the Friedman Eq. (1.3)weare going to split the proof of inequality (3.4) in the following two propositions. Proposition 3.4 There exists λ>0such that λH2(u)≤∞ 0y+1 y P∞(y)(u(x)−u(y))2dxdy:= D(u), (3.9) with u =p/P∞, for all p ∈L1((0,+∞)) ∩L2((0,+∞), P−1 ∞). 123 392 J. A. Cañizo et al. Proof We take 0 <δ<1 and split H2(u)in two parts H2(u)=∞ δ∞ y P∞(x)P∞(y)(u(x)−u(y))2dxdy +δ 0∞ y P∞(x)P∞(y)(u(x)−u(y))2dxdy:= H21(u)+H22(u). For i,j≥0 integers we define Ai,j:= Ii,δ Ij,δ (u(x)−u(y))2dydx=Ii,δ Ij,δ (u(x)−u(y))2dxdy. We can estimate both the left and the right-hand sides of (3.9) by using the quantities Ai,j. Step 1: H21(u)bound.- We start working on the term H21(u)(τ ), where 0 <δ<y< x. By swapping (x,y)in the domain of integration, we get H21(u)=∞ δx δ P∞(x)P∞(y)(u(x)−u(y))2dydx ≤ ∞  i=0 i  j=0Ii,δ Ij,δ (u(x)−u(y))2P∞(x)P∞(y)dydx. Now, using the inequality (3.5) and the symmetry Ai,j=Aj,i, we obtain H21(u)≤B(δ)2 ∞  i=0 i  j=0 pipjIi,δ Ij,δ (u(x)−u(y))2dydx =B(δ)2 ∞  i=0 i  j=0 pipjAi,j =B(δ)2 ∞  j=0 ∞  i=j pipjAi,j=B(δ)2 ∞  i=0 ∞  j=i pipjAi,j.(3.10) Note that some terms in this expression already appear in the right hand side of (3.9), since: ∞  i=0 p2 iAi,i= ∞  i=0 pipiIi,δ Ii,δ (u(x)−u(y))2dxdy ≤1 A(δ)2 ∞  i=0Ii,δ Ii,δ (u(x)−u(y))2P∞(x)P∞(y)dxdy 123 Exponential equilibration of genetic circuits using entropy methods 393 =2 A(δ)2 ∞  i=0Ii,δ x∈Ii,δ x>y (u(x)−u(y))2P∞(x)P∞(y)dxdy ≤2 A(δ)2 ∞  i=0Ii,δ y+1 y (u(x)−u(y))2P∞(x)P∞(y)dxdy =2 A(δ)2∞ δy+1 y (u(x)−u(y))2P∞(x)P∞(y)dxdy ≤PM A(δ)2∞ δy+1 y (u(x)−u(y))2P∞(y)dxdy≤PM A(δ)2D(u), (3.11) where PM=max x∈[δ, ∞)P∞(x)<∞due to the properties described in Sect. 2.1. In order to estimate Ai,jfor j>iwe fix i,jand call n:= j−i≥1. We use n−1 “intermediate reactions” to write the following: introduce n−1 dummy integration variables zi+1,...,zj−1and denote averaged integrals with a stroke. Thus, we have: 4Ai,j=− Ii,δ − Ij,δ (u(x)−u(y))2dxdy =− Ii,δ − Ii+1,δ ···− Ij,δ (u(x)−u(y))2dxdzj−1···dzi+1dy =− Ii,δ − Ii+1,δ ···− Ij,δ u(zj)−u(zi)2dzjdzj−1···dzi, where the last step is just renaming x≡zjand y≡zi. Observe that nothing has been done in the case j=i+1. Using the Cauchy-Schwarz inequality and (3.7), we have 4Ai,j=− Ii,δ − Ii+1,δ ···− Ij,δ ⎛ ⎝ j−1  k=i (u(zk+1)−u(zk))⎞ ⎠ 2 dzjdzj−1···dzi ≤− Ii,δ − Ii+1,δ ···− Ij,δ ⎛ ⎝ j−1  k=i (u(zk+1)−u(zk))2mk⎞ ⎠⎛ ⎝ j−1  k=i 1 mk⎞ ⎠dzjdzj−1···dzi ≤Mj− Ii,δ − Ii+1,δ ···− Ij,δ ⎛ ⎝ j−1  k=i (u(zk+1)−u(zk))2mk⎞ ⎠dzjdzj−1···dzi =Mj j−1  k=i mk− Ii,δ − Ii+1,δ ···− Ij,δ (u(zk+1)−u(zk))2dzjdzj−1···dzi =Mj j−1  k=i mk− Ik,δ − Ik+1,δ (u(zk+1)−u(zk))2dzk+1dzk=4Mj j−1  k=i mkAk,k+1. 123 394 J. A. Cañizo et al. Hence, we deduce that Ai,j≤Mj j−1  k=i mkAk,k+1for all j>i. Thus, we get ∞  i=0 ∞  j=i+1 pipjAi,j≤ ∞  i=0 ∞  j=i+1 pipjMj j−1  k=i mkAk,k+1 = ∞  k=0 mkAk,k+1 ∞  j=k+1 pjMj k  i=0 pi ≤C1 δ ∞  k=0 Ak,k+1mk ∞  j=k+1 Mjpj. The inequality k i=0pi≤C, in the previous expression, holds because ∞ i=0piis a convergent series due to the d’Alembert’s ratio test. Moreover, (3.8) implies ∞  i=0 ∞  j=i+1 pipjAi,j≤C ∞  k=0 Ak,k+1pk(3.12) for a generic constant C>0. We finally work in the Eq. (3.12) to obtain ∞  k=0 Ak,k+1pk= ∞  k=0Ik,δ Ik+1,δ (u(x)−u(y))2dxp kdy ≤1 A(δ) ∞  k=0Ik,δ y+1 y (u(x)−u(y))2dxP ∞(y)dy ≤1 A(δ) ∞ 0y+1 y (u(x)−u(y))2P∞(y)dxdy=1 A(δ) D(u), where we use that y<δ+k+1 2<δ+k+2 2<y+1 and (3.5). We conclude by plugging the above estimate in (3.12), which together with Eqs. (3.10) and (3.11) show that λ1H21(u)≤D(u), (3.13) for some constant λ1>0. Step 2: H22(u)bound.- To prove that there exists λ2>0 such that λ2H22(u)≤∞ 0y+1 y P∞(y)(u(x)−u(y))2dxdy, 123 Exponential equilibration of genetic circuits using entropy methods 395 we use an intermediate variable z∈(δ, 1)as follows: δ 0∞ y (u(x)−u(y))2P∞(x)P∞(y)dxdy =− 1 δδ 0∞ y (u(x)−u(y))2P∞(x)P∞(y)dxdydz ≤2− 1 δδ 0∞ y (u(x)−u(z))2P∞(x)P∞(y)dxdydz +2− 1 δδ 0∞ y (u(z)−u(y))2P∞(x)P∞(y)dxdydz := 2I1+2I2 We bound each of the terms I1,I2.First,forI1we deduce that I1=− 1 δδ 0∞ y (u(x)−u(z))2P∞(x)P∞(y)dxdydz ≤− 1 δ∞ 0 (u(x)−u(z))2P∞(x)dxdz =− 1 δ∞ δ (u(x)−u(z))2P∞(x)dxdz +− 1 δδ 0 (u(x)−u(z))2P∞(x)dxdz:= I11 +I12, since ∞ 0P∞(y)dy=1. For I11 we use that P∞is bounded below on [δ, 1](1 Cδ≤ P∞(x), x∈[δ, 1]) to deduce I11 =− 1 δ∞ δ (u(x)−u(z))2P∞(x)dxdz ≤Cδ− 1 δ∞ δ (u(x)−u(z))2P∞(x)P∞(z)dxdz ≤Cδ 1−δ∞ δ∞ δ (u(x)−u(z))2P∞(x)P∞(z)dxdz =2Cδ 1−δ∞ δ∞ z (u(x)−u(z))2P∞(x)P∞(z)dxdz, Note that the right hand side of the above equation is bounded by a multiple of the term H21(u), thus leading to I11 ≤CH21(u)with C=2Cδ 1−δ.Using(3.13) we deduce that I11 ≤CD(u). 123 396 J. A. Cañizo et al. The integral I12 is clearly smaller than the right hand side of (3.9) since it involves a smaller domain of integration, indeed we obtain I12 =− 1 δδ 0 (u(x)−u(z))2P∞(x)dxdz =− 1 δδ 0 (u(x)−u(z))2P∞(z)dzdx =1 1−δδ 01 δ (u(x)−u(z))2P∞(z)dxdz ≤1 1−δδ 0z+1 z (u(x)−u(z))2P∞(z)dxdz≤CD(u), since z<δ<x<1<z+1. For I2(τ), notice that I2=− 1 δδ 0 (u(z)−u(y))2P∞(y)∞ y P∞(x)dxdydz ≤− 1 δδ 0 (u(z)−u(y))2P∞(y)dydz=I12, and thus, we also deduce that I2≤CD(u). Putting together the estimates on I11,I12 and I2, we conclude that λ2H22(u)≤D(u), (3.14) for some λ2>0. Finally, inequalities (3.13) and (3.14) together imply that λH2(u)≤ D(u)concluding the proof.  Proposition 3.5 There exists α>0such that αD(u)≤D2(u). (3.15) with u =p/P∞, for all p ∈L1((0,+∞)) ∩L2((0,+∞), P−1 ∞). Proof Note that, y<x<y+1 on the left hand side of (3.15). Thus, we can bound the term ω(x−y)with x∈[y,y+1]. Since ω(x)is a decreasing function of x, then ω(1)=1 be−1 b≤ω(x−y)≤1 b=ω(0)with x∈[y,y+1]and y∈R+. Moreover, the term c(x)is bounded, ε≤c(x)≤1 for all x∈R+. So that: ∞ 0y+1 y P∞(y)(u(x)−u(y))2dxdy ≤b εe1 b∞ 0y+1 y ω(x−y)c(y)P∞(y)(u(x)−u(y))2dxdy 123 Exponential equilibration of genetic circuits using entropy methods 397 A 0 20 40 60 80 100 0 0.1 0.2 0.3 0.4 x P∞(x) B 012345 10-2 10-1 100 101 102 Fig. 4 Case 1 Fig 3:H=−4,ε=0.15,K=45,a=5,b=10 ≤b εe1 b∞ 0∞ y ω(x−y)c(y)P∞(y)(u(x)−u(y))2dxdy =b aεe1 bD2(u), which proves the inequality (3.15).  Proof of Theorem 1.1 Putting together (3.9) and (3.15) from Propositions 3.4 and 3.5, we deduce that the entropy-entropy production inequality (3.4) holds. Lemma 3.1 together with (3.4) finally implies (3.2). As consequence, we deduce the exponential convergence towards P∞for all mild solutions of (1.3).  3.2 Numerical illustration of exponential convergence The entropy functional, G2(u)(t), is represented in the plots B of Figs. 4,5,6,7and 8, which address the five possible steady states plots A of Figs. 4,5,6,7and 8(see also Fig. 3). For all cases, these functions are represented in a semi-logarithm scale to numerically validate the exponential convergence shown in the previous section. A Gaussian distribution with mean 2 and standard deviation 0.1, N(2,0.1), has been considered as initial condition. 4ThenD PIDE model We can generalise the entropy functional (2.3) defined for the one dimension PIDE model in order to study the convergence of the multidimensional model. A wellposedness theory of mild and classical solutions satisfying the positivity and mass preservation, the L1-contraction principle, and the maximum principle can be analogously obtained from the one dimensional strategy in Sect. 2. Let us summarize these properties in the next proposition. 123 404 J. A. Cañizo et al. Proof By expanding the square, we can write Gn 2(u)=1 2Rn +Rn + P∞(x)P∞(y)(u(t,x)−u(y))2dxdy.(4.8) We split the latter integral in two parts: the integral over δ×δ, and the integral over its complement with δ= ntimes  ! " [δ, 1/δ]×···×[δ, 1/δ]such that, δ∈(0,1). For the integral over the complement, using p≤C1P∞, we deduce R2n +\(δ×δ) P∞(x)P∞(y)(u(x)−u(y))2dxdy ≤2C2 1R2n +\(δ×δ) P∞(x)P∞(y)dxdy. On the other hand, for the integral over δ×δwe get δδ P∞(x)P∞(y)(u(x)−u(y))2dxdy≤Kδ,1δδ (u(x)−u(y))2dxdy, where Kδ,1:= sup (x,y)∈δ×δ P∞(x)P∞(y)<+∞. We now rewrite u(x)−u(y)as a sum of nterms, each of which being a difference of values of uat points which differ only by one coordinate u(x)−u(y)= n  i=1u(x1,...,xi,yi+1,...,yn)−u(x1,...,xi−1,yi,...,yn), (where it is understood that u(x1,...,xi,yi+1,...,yn)=u(x)for i=n, and u(x1,...,xi−1,yi,...,yn)=u(y)for i=1). Then, by Cauchy-Schwarz’s inequality we have δδ (u(x)−u(y))2dxdy ≤n n  i=1δδu(x1,...,xi,yi+1,...,yn)−u(x1,...,xi−1,yi,...,yn)2 dxdy =n1 δ−δn−1n  i=1[δ,1/δ]n[δ,1/δ]u(x)−u(yi)2 dxidyi 123 Exponential equilibration of genetic circuits using entropy methods 405 =2n1 δ−δn−1n  i=1[δ,1/δ]n1/δ yiu(x)−u(yi)2 dxidyi ≤Kδ,2 n  i=1 ki m[δ,1/δ]n1/δ yi ωi(xi−yi)ci(yi)P∞(yi)u(x)−u(yi)2 dxidyi, therefore we conclude that δδ (u(x)−u(y))2dxdy≤Kδ,2D2(p), (4.9) where Kδ,2is defined by 2n1 δ−δn−1 K−1 δ,2=inf ki mωi(xi−yi)ci(yi)P∞(yi), with the infimum running over all i=1,...,nand over all the points in the domain of integration. We notice that the first of the equalities in (4.9) is just obtained by integrating in the variables that do not appear in the expression and renaming the others; and the second equality is due to the symmetry of the integrand in the variables (xi,yi).Using(4.8)–(4.9) finally gives: Gn 2(u)≤C2 1R2n +\(δ×δ) P∞(x)P∞(y)dxdy+1 2Kδ,1Kδ,2Dn 2(u). We may choose δ>0 such that the first term is smaller than . This gives then the result with K=1 2Kδ,1Kδ,2. Theorem 4.5 (Long-time behaviour) Given any mild solution p with normalised nonnegative initial data p0∈L1(R+)to Eq. (1.8)and given a stationary solution P∞(x) to (1.8)satisfying Assumption 4.1, then lim t→∞ Rn + |p(t,x)−P∞(x)|2dx=0. As a consequence, stationary solutions P∞(x)of (1.8)satisfying Assumption 4.1,if they exist, they are unique. Proof Step 1: Proof for “nice” initial data. We first prove the result for initial data p0∈L1(Rn +)∩C2(Rn +)such that p0≤C1P∞, for some constant C1>0. Observe that this implies in particular that p0∈L2(Rn +,P∞(x)−1dx). For such initial data we deduce that for all t≥0 p(t,x)≤C1P∞(x)for almost all x∈Rn +, 123 406 J. A. Cañizo et al. from the maximum principle. This enables us to use Lemma 4.4. Using the general entropy identity with H(u)=(u−1)2, from Proposition 4.3 we obtain: dGn 2(u) dt=−Dn 2(u). (4.10) Next, by using time integration on [0,T]in Eq. (4.10), the following equality holds for all T>0: Gn 2(u)(T)+T 0 Dn 2(p)(t)dt=Gn 2(u)(0), from which we deduce that: ∞ 0 Dn 2(u)(t)dt<∞.(4.11) From (4.11), there exists a sequence (ts)s≥1such that Dn 2(u)(ts)→0ass→+∞. Thus if we take any >0, then Lemma 4.4 gives: Gn 2(u)(ts)≤KDn 2(u)(ts)+→as s→+∞. Since Gn 2(u)(t)is decreasing in t, this shows that limt→+∞ Gn 2(u)(t)≤. Since is arbitrary chosen, we deduce that: Gn 2(u)(t)→0ast→+∞. Step 2: Proof for all integrable initial data. It is now classical to extend the result in step 1 to all initial data in L1(Rn +)by the L1-contraction principle. In fact, any p0∈L1(Rn +)can be approximated in L1(Rn +)by a sequence (ps 0)s≥1such that ps 0≤sP ∞, for all s≥1. Thus consider the solution psassociated to initial data ps 0. By step 1, we get ∞ 0 |ps(t,x)−P∞(x)|dx→0ast→+∞, since Gn 2(us)(t)≥ps(t,x)−P∞(x)2 1with us=ps P∞. Hence, for s≥1 we deduce ∞ 0 |p(t,x)−P∞(x)|dx≤∞ 0 |p(t,x)−ps(t,x)|dx+∞ 0 |ps(t,x)−P∞(x)|dx ≤∞ 0 |p0(x)−ps 0(x)|dx+∞ 0 |ps(t,x)−P∞(x)|dx, from the L1-contraction principle. This easily leads to the result since lim s→∞ ∞ 0 |p0(x)−ps 0(x)|dx=0 and lim t→∞ ∞ 0 |ps(t,x)−P∞(x)|dx=0, for all s≥1.  123 Exponential equilibration of genetic circuits using entropy methods 407 AB 012345 10-2 10-1 100 101 102 Fig. 9 Example of two self regulated proteins whose distribution has a peak in x=(0,0). (Same parameters as in the example depicted in Fig. 4for both proteins) 4.2 Numerical exploration of the convergence rates The entropy functional, Gn 2(u)(t), is represented in the plots B of Figs. 9,10 and 11, which address three possible steady states (plots A of Figs. 9,10 and 11) that have been obtained using the SELANSI toolboox (Pájaro et al. 2018). For all cases, these functions are represented in a semi-logarithm scale to numerically check if the convergence shown in the previous section is exponential in higher dimensions. In the first example, Fig. 9, we consider two different self-regulated proteins with input functions: ci(xi)=KHi i+εixHi i KHi i+xHi i for i=1,2, with Hi=−4, εi=0.15, Ki=45, ai=5 and bi=10 as in the example depicted in Fig. 4. The second example, Fig. 10, is a self and cross-regulated gene network expressing two different proteins where the first one activates the production of both itself and the second protein, while the second protein inhibits the expression of both proteins. The input functions considered, as in Pájaro et al. (2017), read: c1(x)=11xH11 1xH12 2+12 KH11 11 xH12 2+13xH11 1KH12 12 +KH11 11 KH12 12 xH11 1xH12 2+KH11 11 xH12 2+xH11 1KH12 12 +KH11 11 KH12 12 , c2(x)=21xH22 2xH21 1+22 KH22 22 xH21 1+23xH22 2KH21 21 +KH22 22 KH21 21 xH22 2xH21 1+KH22 22 xH21 1+xH22 2KH21 21 +KH22 22 KH21 21 , (4.12) with H11 =−4, H21 =−6, H12 =H22 =2, K11 =K12 =45, K21 =K22 =70, ε11 =ε21 =0.002, ε12 =0.02, ε22 =0.1, ε13 =ε23 =0.2 and network parameters γ1 x=γ2 x=1, γ1 m=γ2 m=25, k1 m=10, k2 m=20, b1=10 and b2=20. 123 408 J. A. Cañizo et al. AB 012345 10-5 100 Fig. 10 Example of two self and cross regulated proteins whose distribution has a peak in some positive point x=(x1,x2)with x1>0andx2>0. Parameters: γ1 x=γ2 x=1, γ1 m=γ2 m=25, k1 m=10, k2 m=20, b1=10, b2=20 and input functions in (4.12) AB 0246810 10-10 10-5 100 105 1010 Fig. 11 Example of two mutual repressed proteins whose joint distribution is bimodal attaining two peaks in two positive points. Parameters:γ1 x=γ2 x=1, γ1 m=γ2 m=25, k1 m=k2 m=8andb1=b2=16 with input functions defined in (4.13) Our third example, Fig. 11, corresponds to a mutual repressing network of two genes in which the protein produced by the expression of one gene inhibits the production of the other protein in the network. The input functions, as in Pájaro et al. (2017), for this example take the following form: c1(x)=KH12 1+ε1xH12 2 KH12 1+xH12 2 ,c2(x)=KH21 2+ε2xH21 1 KH21 2+xH21 1 ,(4.13) with H12 =H21 =4, K1=K2=45 and ε1=ε2=0.15. The dimensionless network parameters are γ1 x=γ2 x=1, γ1 m=γ2 m=25, k1 m=k2 m=8 and b1=b2= 16. 123 Exponential equilibration of genetic circuits using entropy methods 409 For each example described above, a multivariate Gaussian distribution with means 10 and standard deviations 1, N([10,10],[1,1]), has been considered as initial condition. 5 Conclusions Analytical results for the nDmodel show convergence to equilibrium via a very general method, but do not give a bound on the convergence rate. The numerical simulations we have carried out clearly support the idea that exponential convergence also holds in the multidimensional case, though we have not been able to prove this using the same entropy method as in the one-dimensional case. Approach to equilibrium seems to follow a steady exponential speed, being quickly dominated by the spectral gap expected from our analysis. There also seem to be initial regimes where the approach to equilibrium can occur much faster; our interpretation is that smaller (more negative) eigenvalues can dominate at initial stages of time evolution, but are overcome by the dominant eigenvalue as equilibrium is approached. Acknowledgements J. A. Cañizo and J. A. Carrillo were supported by Projects MTM2014-52056-P and MTM2017-85067-P, funded by the Spanish government and the European Regional Development Fund. J. A. Carrillo was partially supported by the EPSRC Grant Number EP/P031587/1. M. Pájaro acknowledges support from Spanish MINECO fellowships BES-2013-063112, EEBB-I-16-10540 and EEBB-I-17-12182. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. References Alon U (2007) An introduction to systems biology. Design principles of biological circuits. Chapman & Hall/ CRC, London Balagué D, Cañizo JA, Gabriel P (2013) Fine asymptotics of profiles and relaxation to equilibrium for growth-fragmentation equations with variable drift rates. Kinet Relat Models 6(2):219–243 Bokes P, Singh A (2015) Protein copy number distributions for a self-regulating gene in the presence of decoy binding sites. PLoS ONE 10(3):e0120555 Bokes P, Singh A (2017) Gene expression noise is affected differentially by feedback in burst frequency and burst size. J Math Biol 74(6):1483–1509 Bokes P, Lin YT, Singh A (2018) High cooperativity in negative feedback can amplify noisy gene expression. Bull Math Biol 80(7):1871–1899 Cáceres MJ, Cañizo JA, Mischler S (2011) Rate of convergence to an asymptotic profile for the self-similar fragmentation and growth-fragmentation equations. J Math Pures Appl 96(4):334–362 Cañizo JA, Carrillo JA, Cuadrado SL (2013) Measure solutions for some models in population dynamics. Acta Appl Math 123:141–156 Carrillo JA, Cordier S, Mancini S (2011) A decision-making fokker-planck model in computational neuroscience. J Math Biol 63(5):801–830 Dar RD, Razooky BS, Singh A, Trimeloni TV, McCollum JM, Cox CD, Simpson ML, Weinberger LS (2012) Transcriptional burst frequency and burst size are equally modulated across the human genome. Proc Natl Acad Sci USA 109(43):17454–17459 Doumic Jauffret M, Gabriel P (2010) Eigenelements of a general aggregation-fragmentation model. Math Models Methods Appl Sci 20(5):757–783 123 410 J. A. Cañizo et al. Elgart V, Jia T, Fenley AT, Kulkarni R (2011) Connecting protein and mRNA burst distributions for stochastic models of gene expression. Phys Biol 8:046001 Elowitz MB, Levine AJ, Siggia ED, Swain PS (2002) Stochastic gene expression in a single cell. Science 297(5584):1183–1186 Engblom S (2006) Computing the moments of high dimensional solutions of the master equation. Appl Math Comput 180(2):498–515 Engel K-J, Nagel R (2006) A short course on operator semigroups. Universitext. Springer, New York Friedman N, Cai L, Xie XS (2006) Linking stochastic dynamics to population distribution: an analytical framework of gene expression. Phys Rev Lett 97(16):168302 Gillespie DT (1976) A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys 22(4):403–434 Gillespie DT (2007) Stochastic simulation of chemical kinetics. Annu Rev Phys Chem 58:35–55 Gualdani MP, Mischler S, Mouhot C (2010) Factorization for non-symmetric operators and exponential h-theorem. June Hasenauer J, Wolf V, Kazeroonian A, Theis FJ (2015) Method of conditional moments (mcm) for the chemical master equation: a unified framework for the method of moments and hybrid stochasticdeterministic models. J Math Biol. 69(3):687–735 Jahnke T (2011) On reduced models for the chemical master equation. Multiscale Model. Simul. 9(4):1646– 1676 Kærn M, Elston TC, Blake WJ, Collins JJ (2005) Stochasticity in gene expression: from theories to phenotypes. Nat Rev Genet 6(6):451–464 Kepler TB, Elston TC (2001) Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations. Biophys J 81(6):3116–3136 Laurençot P, Perthame B (2009) Exponential decay for the growth-fragmentation/cell-division equation. Commun Math Sci 7(2):503–510 Lee TH, Maheshri N (2012) A regulatory role for repeated decoy transcription factor binding sites in target gene expression. Mol Syst Biol 8:576 Mackey MC, Tyran-Kaminska M, Yvinec R (2011) Molecular distributions in gene regulatory dynamics. J Theor Biol 274(1):84–96 McAdams H, Arkin A (1997) Stochastic mechanisms in gene expression. Proc Natl Acad Sci USA 94:814– 819 Michel P, Mischler S, Perthame B (2004) General entropy equations for structured population models and scattering. C R Math 338(9):697–702 Michel P, Mischler S, Perthame B (2005) General relative entropy inequality: an illustration on growth models. J Math Pures Appl 84(9):1235–1260 Munsky B, Khammash M (2006) The finite state projection algorithm for the solution of the chemical master equation. J Chem Phys 124(4):1–12 Ochab-Marcinek A, Tabaka M (2015) Transcriptional leakage versus noise: a simple mechanism of conversion between binary and graded response in autoregulated genes. Phys Rev E 91(1):012704 Pájaro M, Alonso AA, Vázquez C (2015) Shaping protein distributions in stochastic self-regulated gene expression networks. Phys Rev E 92(3):032712 Pájaro M, Alonso AA, Carrillo JA, Vázquez C (2016) Stability of stochastic gene regulatory networks using entropy methods. IFAC-PapersOnLine 49(24):1–5 Pájaro M, Alonso AA, Otero-Muras I, Vázquez C (2017) Stochastic modeling and numerical simulation of gene regulatory networks with protein bursting. J Theor Biol 421:51–70 Pájaro M, Otero-Muras I, Vázquez C, Alonso AA (2018) SELANSI: a toolbox for simulation of stochastic gene regulatory networks. Bioinformatics 34(5):893–895 Paulsson J (2004) Summing up the noise in gene networks. Nature 427:415–418 Paulsson J (2005) Models of stochastic gene expression. Phys Life Rev 2(2):157–175 Perthame B (2007) Transport equations in biology. Frontiers in mathematics. Birkhäuser Verlag, Basel Perthame B, Ryzhik L (2005) Exponential decay for the fragmentation or cell-division equation. J Differ Equ 210(1):155–177 Shahrezaei V, Swain PS (2008) Analytical distributions for stochastic gene expressions. Proc Natl Acad Sci USA 105(45):17256–17261 Sherman MS, Cohen BA (2014) A computational framework for analyzing stochasticity in gene expression. PLoS Comput Biol 10(5):1003596 123 Exponential equilibration of genetic circuits using entropy methods 411 Thomas P, Popovic N, Grima R (2014) Phenotypic switching in gene regulatory networks. Proc Natl Acad Sci USA 111(19):6994–6999 Van Kampen NG (2007) Stochastic processes in physics and chemistry, 3rd edn. Elsevier, Amsterdam Wallace EWJ, Gillespie DT, Sanft KR, Petzold LR (2012) Linear noise approximation is valid over limited times for any chemical system that is sufficiently large. IET Syst Biol 6(4):102–115 123