Michael C Mackey
email: mackey@cnd.mcgill.ca

Marta Tyran-Kamińska

Romain Yvinec
email: yvinec@math.univ-lyon1.fr

Molecular Distributions in Gene Regulatory Dynamics

Keywords:

We show how one may analytically compute the stationary density of the distribution of molecular constituents in populations of cells in the presence of noise arising from either bursting transcription or translation, or noise in degradation rates arising from low numbers of molecules. We have compared our results with an analysis of the same model systems (either inducible or repressible operons) in the absence of any stochastic effects, and shown the correspondence between behaviour in the deterministic system and the stochastic analogs. We have identified key dimensionless parameters that control the appearance of one or two steady states in the deterministic case, or unimodal and bimodal densities in the stochastic systems, and detailed the analytic requirements for the occurrence of different behaviours. This approach provides, in some situations, an alternative to computationally intensive stochastic simulations. Our results indicate that, within the context of the simple models we have examined, bursting and degradation noise cannot be distinguished analytically when present alone.

Introduction

In neurobiology, when it became clear that some of the fluctuations seen in whole nerve recording, and later in single cell recordings, were not simply measurement noise but actual fluctuations in the system being studied, researchers very quickly started wondering to what extent these fluctuations actually played a role in the operation of the nervous system.

Much the same pattern of development has occurred in cellular and molecular biology as experimental techniques have allowed investigators to probe temporal behaviour at ever finer levels, even to the level of individual molecules. Experimentalists and theoreticians alike who are interested in the regulation of gene networks are increasingly focussed on trying to access the role of various types of fluctuations on the operation and fidelity of both simple and complex gene regulatory systems. Recent reviews [START_REF] Kaern | Stochasticity in gene expression: From theories to phenotypes[END_REF][START_REF] Raj | Nature, nurture, or chance: Stochastic gene expression and its consequences[END_REF]Shahrezaei and Swain, 2008b) give an interesting perspective on some of the issues confronting both experimentalists and modelers.

Typically, the discussion seems to focus on whether fluctuations can be considered as extrinsic to the system under consideration (Shahrezaei et al., 2008;[START_REF] Ochab-Marcinek | Predicting the asymmetric response of a genetic switch to noise[END_REF][START_REF] Ochab-Marcinek | Extrinsic noise passing through a Michaelis-Menten reaction: A universal response of a genetic switch[END_REF], or whether they are an intrinsic part of the fundamental processes they are affecting (e.g. bursting, see below). The dichotomy is rarely so sharp however, but [START_REF] Elowitz | Stochastic gene expression in a single cell[END_REF] have used an elegant experimental technique to distinguish between the two, see also [START_REF] Raser | Control of stochasticity in eukaryotic gene expression[END_REF], while [START_REF] Swain | Intrinsic and extrinsic contributions to stochasticity in gene expression[END_REF] and [START_REF] Scott | Estimations of intrinsic and extrinsic noise in models of nonlinear genetic networks[END_REF] have laid the groundwork for a theoretical consideration of this question. One issue that is raised persistently in considerations of the role of fluctuations or noise in the operation of gene regulatory networks is whether or not they are "beneficial" [START_REF] Blake | Phenotypic consequences of promotermediated transcriptional noise[END_REF] or "detrimental" [START_REF] Fraser | Noise minimization in eukaryotic gene expression[END_REF] to the operation of the system under consideration. This is, of course, a question of definition and not one that we will be further concerned with here.

Here, we consider in detail the density of the molecular distributions in generic bacterial operons in the presence of 'bursting' (commonly known as intrinsic noise in the biological literature) as well as inherent (extrinsic) noise using an analytical approach. Our work is motivated by the well documented production of mRNA and/or protein in stochastic bursts in both prokaryotes and eukaryotes [START_REF] Blake | Noise in eukaryotic gene expression[END_REF][START_REF] Cai | Stochastic protein expression in individual cells at the single molecule level[END_REF][START_REF] Chubb | Transcriptional pulsing of a developmental gene[END_REF][START_REF] Golding | Real-time kinetics of gene activity in individual bacteria[END_REF][START_REF] Raj | Stochastic mRNA synthesis in mammalian cells[END_REF][START_REF] Sigal | Variability and memory of protein levels in human cells[END_REF][START_REF] Yu | Probing gene expression in live cells, one protein molecule at a time[END_REF], and follows other contributions by, for example, [START_REF] Kepler | Stochasticity in transcriptional regulation: Origins, consequences, and mathematical representations[END_REF], [START_REF] Friedman | Linking stochastic dynamics to population distribution: An analytical framework of gene expression[END_REF], [START_REF] Bobrowski | Asymptotic behavior of distributions of mRNA and protein levels in a model of stochastic gene expression[END_REF] and Shahrezaei and Swain (2008a).

In Section 2 we develop the concept of the operon and treat simple models of the classic inducible and repressible operon. Section 4 considers the effects of bursting alone in an ensemble of single cells. Section 5 then examines the situation in which there are continuous white noise fluctuations in the dominant species degradation rate in the absence of bursting.

Generic operons

The operon concept

The so-called 'central dogma' of molecular biology is simple to state in principle, but complicated in its detail. Namely through the process of transcription of DNA, messenger RNA (mRNA, M) is produced and, in turn, through the process of translation of the mRNA proteins (intermediates, I) are produced. There is often feedback in the sense that molecules (enzymes, E) whose production is controlled by these proteins can modulate the translation and/or transcription processes. In what follows we will refer to these molecules as effectors. We now consider both the transcription and translation process in more detail.

In the transcription process an amino acid sequence in the DNA is copied by the enzyme RNA polymerase (RNAP) to produce a complementary copy of the DNA segment encoded in the resulting RNA. Thus this is the first step in the transfer of the information encoded in the DNA. The process by which this occurs is as follows.

When the DNA is in a double stranded configuration, the RNAP is able to recognize and bind to the promoter region of the DNA. (The RNAP/double stranded DNA complex is known as the closed complex.) Through the action of the RNAP, the DNA is unwound in the vicinity of the RNAP/DNA promoter site, and becomes single stranded. (The RNAP/single stranded DNA is called the open complex.) Once in the single stranded configuration, the transcription of the DNA into mRNA commences.

In prokaryotes, translation of the newly formed mRNA commences with the binding of a ribosome to the mRNA. The function of the ribosome is to 'read' the mRNA in triplets of nucleotide sequences (codons). Then through a complex sequence of events, initiation and elongation factors bring transfer RNA (tRNA) into contact with the ribosome-mRNA complex to match the codon in the mRNA to the anti-codon in the tRNA. The elongating peptide chain consists of these linked amino acids, and it starts folding into its final conformation. This folding continues until the process is complete and the polypeptide chain that results is the mature protein.

The lactose (lac) operon in bacteria is the paradigmatic example of this concept and this much studied system consists of three structural genes named lacZ, lacY, and lacA. These three genes contain the code for the ultimate production, through the translation of mRNA, of the intermediates β-galactosidase, lac permease, and thiogalactoside transacetylase respectively. The enzyme β-galactosidase is active in the conversion of lactose into allolactose and then the conversion of allolactose into glucose. The lac permease is a membrane protein responsible for the transport of extracellular lactose to the interior of the cell. (Only the transacetylase plays no apparent role in the regulation of this system.) The regulatory gene lacI, which is part of a different operon, codes for the lac repressor, which is transformed to an inactive form when bound with allolactose, so in this system allolactose functions as the effector molecule.

The transcription rate function

In this section we examine the molecular dynamics of both the classical inducible and repressible operon to derive expressions for the dependence of the transcription rate on effector levels. (When the transcription rate is constant and independent of the effector levels we will refer to this as the no control situation.)

Inducible regulation

For a typical inducible regulatory situation (such as the lac operon), in the presence of the effector molecule the repressor is inactive (is unable to bind to the operator region preceding the structural genes), and thus DNA transcription can proceed. Let R denote the repressor, E the effector molecule, and O the operator. The effector is known to bind with the active form R of the repressor. We assume that this reaction is of the form

R + nE K 1 ⇋ RE n K 1 = RE n R • E n , ( 1
)
where n is the effective number of molecules of effector required to inactivate the repressor R. Furthermore, the operator O and repressor R are assumed to interact according to

O + R K 2 ⇋ OR K 2 = OR O • R .
Let the total operator be O tot :

O tot = O + OR = O + K 2 O • R = O(1 + K 2 R),
and the total level of repressor be R tot :

R tot = R + K 1 R • E n + K 2 O • R.
The fraction of operators not bound by repressor (and therefore free to synthesize mRNA) is given by

f (E) = O O tot = 1 1 + K 2 R .
If the amount of repressor R bound to the operator O is small

R tot ≃ R + K 1 R • E n = R(1 + K 1 E n ) so R = R tot 1 + K 1 E n , 2 and consequently f (E) = 1 + K 1 E n 1 + K 2 R tot + K 1 E n = 1 + K 1 E n K + K 1 E n , (2)
where K = 1 + K 2 R tot . There will be maximal repression when E = 0 but even then there will still be a basal level of mRNA production proportional to K -1 (which we call the fractional leakage).

If the maximal DNA transcription rate is φm (in units of inverse time) then, under the assumption that the rate of transcription ϕ in the entire population is proportional to the fraction f of unbound operators, the variation ϕ of the DNA transcription rate with the effector level is given by ϕ = φm f , or

ϕ(E) = φm 1 + K 1 E n K + K 1 E n .
(3)

Repressible regulation

In the classic example of a repressible system (such as the trp operon) in the presence of the effector molecule the repressor is active (able to bind to the operator region), and thus block DNA transcription. We use the same notation as before, but now note that the effector binds with the inactive form R of the repressor so it becomes active. We assume that this reaction is of the same form as in Equation 1. The difference now is that the operator O and repressor R are assumed to interact according to

O + R • E n K 2 ⇋ ORE n K 2 = ORE n O • R • E n .
The total operator is now given by

O tot = O + ORE n = O + K 2 O • R • E n = O(1 + K 2 R • E n ),
so the fraction of operators not bound by repressor is given by

f (E) = O O tot = 1 1 + K 2 R • E n .
Again assuming that the amount of repressor R bound to the operator O is small we have

f (E) = 1 + K 1 E n 1 + (K 1 + K 2 R tot )E n = 1 + K 1 E n 1 + KE n ,
where K = K 1 + K 2 R tot . Now there will be maximal repression when E is large, but even at maximal repression there will still be a basal level of mRNA production proportional to K 1 K -1 < 1. The variation of the DNA transcription rate with effector level is given by ϕ = φm f or

ϕ(E) = φm 1 + K 1 E n 1 + KE n . (4)
Both ( 3) and ( 4) are special cases of the function

ϕ(E) = φm 1 + K 1 E n A + BE n = φm f (E). ( 5
)
where A, B ≥ 0 are given in Table 1. parameter inducible repressible

A K = 1 + K 2 R tot 1 B K 1 K = K 1 + K 2 R tot B A K 1 K K Λ = A K 1 ∆ = BK -1 1 1 KK -1 1 θ = κ d n∆ 1 - ∆ Λ κ d n • K -1 K > 0 κ d n • K 1 -K K < 0

Deterministic operon dynamics in a population of cells

The reader may wish to consult [START_REF] Polynikis | Comparing differeent ODE modelling approaches for gene regulatory networks[END_REF] for an interesting survey of techniques applicable to this approach.

We first consider a large population of cells, each of which contains one copy of a particular operon, and let (M, I, E) denote mRNA, intermediate protein, and effector levels respectively in the population. Then for a generic operon with a maximal level of transcription bd (in concentration units), we have dynamics described by the system (Griffith, 1968a,b;[START_REF] Othmer | The qualitative dynamics of a class of biochemical control circuits[END_REF][START_REF] Selgrade | Mathematical analysis of a cellular control process with positive feedback[END_REF]

dM dt = bd φm f (E) -γ M M, ( 6
)
dI dt = β I M -γ I I, ( 7
)
dE dt = β E I -γ E E. ( 8
)
Here we assume that the rate of mRNA production is proportional to the fraction of time the operator region is active, and that the rates of intermediate and enzyme production are simply proportional to the amount of mRNA and intermediate respectively. All three of the components (M, I, E) are subject to random loss. The function f is calculated in the previous section.

It will greatly simplify matters to rewrite Equations 6-8 by defining dimensionless concentrations. To this end we first rewrite Equation 5 in the form

ϕ(e) = ϕ m f (e), (9)
where ϕ m (dimensionless) is defined by

ϕ m = φm β E β I γ M γ E γ I and f (e) = 1 + e n Λ + ∆e n , (10)
Λ and ∆ are defined in Table 1, and we have defined a dimensionless effector concentration (e) through

E = ηe with η = 1 n √ K 1 .
Further defining dimensionless intermediate (i) and mRNA concentrations (m) through

I = iη γ E β E and M = mη γ E γ I β E β I ,
Equations 6-8 can be written in the equivalent form

dm dt = γ M [κ d f (e) -m], di dt = γ I (m -i), de dt = γ E (i -e),
where

κ d = b d ϕ m and b d = bd η (11)
are dimensionless constants.

For notational simplicity, henceforth we denote dimensionless concentrations by (m, i, e) = (x 1 , x 2 , x 3 ), and subscripts (M, I, E) = (1, 2, 3). Thus we have

dx 1 dt = γ 1 [κ d f (x 3 ) -x 1 ], ( 12
)
dx 2 dt = γ 2 (x 1 -x 2 ), ( 13
)
dx 3 dt = γ 3 (x 2 -x 3 ). (14)
In each equation, γ i for i = 1, 2, 3 denotes a net loss rate (units of inverse time), and thus Equations 12-14 are not in dimensionless form. The dynamics of this classic operon model can be fully analyzed. Let X = (x 1 , x 2 , x 3 ) and denote by S t (X) the flow generated by the system (12)-( 14). For both inducible and repressible operons, for all initial conditions X 0 = (x 0 1 , x 0 2 , x 0 3 ) ∈ R + 3 the flow S t (X 0 ) ∈ R + 3 for t > 0. Steady states of the system (12)-( 14) are in a one to one correspondence with solutions of the equation

x κ d = f (x) (15)
and for each solution x * of Equation 15 there is a steady state

X * = (x * 1 , x * 2 , x *
3 ) of ( 12)-( 14) given by

x * 1 = x * 2 = x * 3 = x * .
Whether there is a single steady state X * or there are multiple steady states will depend on whether we are considering a repressible or inducible operon.

No control

In this case, f (x) ≡ 1, and there is a single steady state x * = κ d that is globally asymptotically stable. 10, and the straight lines correspond to x/κ d for (in a clockwise direction)

κ d ∈ [0, κ d-), κ d = κ d-, κ d ∈ (κ d-, κ d+ ), κ d = κ d+ ,
and κ d+ < κ d . This figure was constructed with n = 4 and K = 10 for which κ d-= 3.01 and κ d+ = 5.91 as computed from (18). See the text for further details.

Inducible regulation

Single versus multiple steady states. For an inducible operon with f given by Equation 2, there may be one (X * 1 or

X * 3 ), two (X * 1 , X * 2 = X * 3 or X * 1 = X * 2 , X * 3 ), or three (X * 1 , X * 2 , X *
3 ) steady states, with the ordering 0 < X * 1 ≤ X * 2 ≤ X * 3 , corresponding to the possible solutions of Equation 15 (cf. Figure 1). The smaller steady state (X * 1 ) is typically referred to as an uninduced state, while the largest steady state (X *

3 ) is called the induced state. The steady state values of x are easily obtained from (15) for given parameter values, and the dependence on κ d for n = 4 and a variety of values of K is shown in Figure 1. Figure 2 shows a graph of the steady states x * versus κ d for various values of the leakage parameter K.

Analytic conditions for the existence of one or more steady states can be obtained by using Equation 15 in conjunction with the observation that the delineation points are marked by the values of κ d at which x/κ d is tangent to f (x) (see Figure 1). Simple differentiation of (15) yields the second condition

1 κ d n(K -1) = x n-1 (K + x n ) 2 . ( 16
)
From equations ( 15) and ( 16) we obtain the values of x at which tangency will occur:

x ± = n K -1 2        n - K + 1 K -1 ± n 2 -2n K + 1 K -1 + 1        . (17)
The two corresponding values of κ d at which a tangency occurs are given by

κ d± = x ∓ K + x n ∓ 1 + x n ∓ . (18)
(Note the deliberate use of x ∓ as opposed to x ± .) A necessary condition for the existence of two or more steady states is obtained by requiring that the square root in (17) be non-negative, or

K ≥ n + 1 n -1 2 . ( 19
)
From this a second necessary condition follows, namely

κ d ≥ n + 1 n -1 n n + 1 n -1 . ( 20
)
Further, from Equations 15 and 16 we can delineate the boundaries in (K, κ d ) space in which there are one or three locally stable steady states as shown in Figure 3. There, we have given a parametric plot (x is the parameter) of κ d versus K, using

K(x) = x n [x n + (n + 1)] (n -1)x n -1 and κ d (x) = [K(x) + x n ] 2 nx n-1 [K(x) -1] ,
for n = 4 obtained from Equations 15 and 16. As is clear from the figure, when leakage is appreciable (small K, e.g for n = 4, K < (5/3) 2 ) then the possibility of bistable behaviour is lost.

Remark 1. Some general observations on the influence of n, K, and κ d on the appearance of bistability in the deterministic case are in order.

1. The degree of cooperativity (n) in the binding of effector to the repressor plays a significant role. Indeed, n > 1 is a necessary condition for bistability. 2. If n > 1 then a second necessary condition for bistability is that K satisfies Equation 19 so the fractional leakage (K -1 ) is sufficiently small. 3. Furthermore, κ d must satisfy Equation 20which is quite instructive. Namely for n → ∞ the limiting lower limit is κ d > 1 while for n → 1 the minimal value of κ d becomes quite large. This simply tells us that the ratio of the product 1 , X * 3 , while outside there is only one. The tip of the cone occurs at (K, κ d ) = ((5/3) 2 , (5/3) 4 √ 5/3) as given by Equations 19 and 20. For K ∈ [0, (5/3) 2 ) there is but a single steady state. of the production rates to the product of the degradation rates must always be greater than 1 for bistability to occur, and the lower the degree of cooperativity (n) the larger the ratio must be. 4. If n, K and κ d satisfy these necessary conditions then bistability is only possible if κ d ∈ [κ d-, κ d+ ] (c.f. Figure 3). 5. The locations of the minimal (x -) and maximal (x + ) values of x bounding the bistable region are independent of κ d .

6. Finally (a) (x + -x -) is a decreasing function of increasing n for constant κ d , K (b) (x + -x -) is an increasing function of increasing K for constant n, κ d .
Local and global stability. The local stability of a steady state x * is determined by the solutions of the eigenvalue equation [START_REF] Yildirim | Dynamics and bistability in a reduced model of the lac operon[END_REF]) 21) Set

(λ + γ 1 )(λ + γ 2 )(λ + γ 3 ) -γ 1 γ 2 γ 3 κ d f ′ * = 0, f ′ * = f ′ (x * ). (
a 1 = 3 i=1 γ i , a 2 = 3 i j=1 γ i γ j , a 3 = (1 -κ d f ′ * ) 3 i=1 γ i , so (21
) can be written as

λ 3 + a 1 λ 2 + a 2 λ + a 3 = 0. ( 22
)
By Descartes's rule of signs, (22) will have either no positive roots for f ′ * ∈ [0, κ -1 d ) or one positive root otherwise. With this information and using the notation SN to denote a locally stable node, HS a half or neutrally stable steady state, and US an unstable steady state (saddle point), then there will be: • A single steady state

X * 1 (SN), for κ d ∈ [0, κ d-)
• Two coexisting steady states X * 1 (SN) and X * 2 = X * 3 (HS, born through a saddle node bifurcation) for

κ d = κ d- • Three coexisting steady states X * 1 (S N), X * 2 (US ), X * 3 (SN) for κ d ∈ (κ d-, κ d+ )
• Two coexisting steady states X * 1 = X * 2 (HS at a saddle node bifurcation), and X * 3 (SN) for κ d = κ d+

• One steady state X * 3 (SN) for κ d+ < κ d .

For the inducible operon, other work extends these local stability considerations and we have the following result characterizing the global behaviour: Theorem 1. [START_REF] Othmer | The qualitative dynamics of a class of biochemical control circuits[END_REF]Smith, 1995, Proposition 2.1, Chapter 4) For an inducible operon with ϕ given by Equation 3, define

I I = [1/K, 1]. There is an attracting box B I ⊂ R + 3 defined by B I = {(x 1 , x 2 , x 3 ) : x i ∈ I I , i = 1, 2, 3}
such that the flow S t is directed inward everywhere on the surface of B I . Furthermore, all X * ∈ B I and

1. If there is a single steady state, i.e. X * 1 for κ d ∈ [0, κ d-), or X *
3 for κ d+ < κ d , then it is globally stable. 2. If there are two locally stable nodes, i.e. X * 1 and X * 3 for κ d ∈ (κ d-, κ d+ ), then all flows S (X 0 ) are attracted to one of them. (See [START_REF] Selgrade | Mathematical analysis of a cellular control process with positive feedback[END_REF] for a delineation of the basin of attraction of X * 1 and X * 3 .)

Repressible regulation

As illustrated in Figure 4, the repressible operon has a single steady state corresponding to the unique solution x * of Equation 15. To determine its local stability we apply the Routh-Hurwitz criterion to the eigenvalue equation ( 22). The steady state corresponding to x * will be locally stable (i.e. have eigenvalues with negative real parts) if and only if a 1 > 0 (always the case) and a 1 a 2 -a 3 > 0.

(23)

The well known relation between the arithmetic and geometric means 1

n n i=1 γ i ≥        n i=1 γi        1/n
, when applied to both a 1 and a 2 gives, in conjunction with Equation 23,

a 1 a 2 -a 3 ≥ (8 + κ d f ′ * ) 3 i=1 γ i > 0.
Thus as long as f ′ * > -8/κ d , the steady state corresponding to x * will be locally stable. Once condition ( 23) is violated, stability of x * is lost via a supercritical Hopf bifurcation and a limit cycle is born. One may even compute the Hopf period of this limit cycle by assuming that λ = jω H ( j = √ -1) in Equation 22where ω H is the Hopf angular frequency. Equating real and imaginary parts of the resultant yields ω H = √ a 3 /a 1 or

T H = 2π ω H = 2π × 3 i=1 γ i (1 -κ d f ′ * ) 3 i=1 γ i .
These local stability results tell us nothing about the global behaviour when stability is lost, but it is possible to characterize the global behaviour of a repressible operon with the following Theorem 2. (Smith, 1995, Theorem 4.1 & Theorem 4.2, Chapter 3) For a repressible operon with ϕ given by Equation 4, define

I R = [K 1 /K, 1]. There is a globally attracting box B R ⊂ R + 3 defined by B R = {(x 1 , x 2 , x 3 ) : x i ∈ I R , i = 1, 2, 3}
such that the flow S is directed inward everywhere on the surface of B R . Furthermore there is a single steady state X * ∈ B R . If X * is locally stable it is globally stable, but if X * is unstable then a generalization of the Poincare-Bendixson theorem (Smith, 1995, Chapter 3) implies the existence of a globally stable limit cycle in B R .

Remark 2. There is no necessary connection between the Hopf period computed from the local stability analysis and the period of the globally stable limit cycle.

Fast and slow variables

In dynamical systems, considerable simplification and insight into the behaviour can be obtained by identifying fast and slow variables. This technique is especially useful when one is initially interested in the approach to a steady state. In this context a fast variable is one that relaxes much more rapidly to an equilibrium than a slow variable [START_REF] Haken | Synergetics: An introduction[END_REF]. In many systems, including chemical and biochemical ones, this is often a consequence of differences in degradation rates, with the fastest variable the one that has the largest degradation rate. We employ the same strategy here to obtain approximations to the population level dynamics that will be used in the next section.

It is often the case that the degradation rate of mRNA is much greater than the corresponding degradation rates for both the intermediate protein and the effector (γ 1 ≫ γ 2 , γ 3 ) so in this case the mRNA dynamics are fast and we have the approximate relationship

x 1 ≃ κ d f (x 3 ).

Consequently the three variable system describing the generic operon reduces to a two variable one involving the slower intermediate and effector:

dx 2 dt = γ 2 [κ d f (x 3 ) -x 2 ], ( 24
)
dx 3 dt = γ 3 (x 2 -x 3 ). ( 25
)
In our considerations of specific single operon dynamics below we will also have occasion to examine two further subcases, namely Case 1. Intermediate (protein) dominated dynamics. If it should happen that γ 1 ≫ γ 3 ≫ γ 2 (as for the lac operon, then the effector also qualifies as a fast variable so

x 3 ≃ x 2
and thus from ( 24)-( 25) we recover the one dimensional equation for the slowest variable, the intermediate:

dx 2 dt = γ 2 [κ d f (x 2 ) -x 2 ]. ( 26
)
Case 2. Effector (enzyme) dominated dynamics. Alternately, if γ 1 ≫ γ 2 ≫ γ 3 then the intermediate is a fast variable relative to the effector and we have

x 2 ≃ x 3 so our two variable system ( 24)-( 25) reduces to a one dimensional system

dx 3 dt = γ 3 [κ d f (x 3 ) -x 3 ] (27)
for the relatively slow effector dynamics. Both Equations 26 and 27 are of the form

dx dt = γ[κ d f (x) -x] ( 28
)
where γ is either γ 2 for protein (x 2 ) dominated dynamics or γ 3 for effector (x 3 ) dominated dynamics.

Distributions with intrinsic bursting

Generalities

It is well documented experimentally [START_REF] Cai | Stochastic protein expression in individual cells at the single molecule level[END_REF][START_REF] Chubb | Transcriptional pulsing of a developmental gene[END_REF][START_REF] Golding | Real-time kinetics of gene activity in individual bacteria[END_REF][START_REF] Raj | Stochastic mRNA synthesis in mammalian cells[END_REF][START_REF] Sigal | Variability and memory of protein levels in human cells[END_REF][START_REF] Yu | Probing gene expression in live cells, one protein molecule at a time[END_REF] that in some organisms the amplitude of protein production through bursting translation of mRNA is exponentially distributed at the single cell level with density

h(y) = 1 b e -y/ b, ( 29
)
where b is the average burst size, and that the frequency of bursting ϕ is dependent on the level of the effector. Writing Equation 29 in terms of our dimensionless variables we have

h(x) = 1 b e -x/b . ( 30
)
Remark 3. The technique of eliminating fast variables described in Section 2.3 above (also known as the adiabatic elimination technique [START_REF] Haken | Synergetics: An introduction[END_REF]) has been extended to stochastically perturbed systems when the perturbation is a Gaussian distributed white noise, c.f. Stratonovich (1963, Chapter 4, Section 11.1), [START_REF] Wilemski | On the derivation of Smoluchowski equations with corrections in the classical theory of Brownian motion[END_REF], [START_REF] Titular | A systematic solution procedure for the Fokker-Planck equation of a Brownian particle in the high-friction case[END_REF], and Gardiner (1983, Section 6.4). However, to the best of our knowledge, this type of approximation has never been extended to the situation dealt with here in which the perturbation is a jump Markov process.

The single cell analog of the population level intermediate protein dominated Case 1 above (when

γ 1 ≫ γ 3 ≫ γ 2 ) is dx 2 dt = -γ 2 x 2 + Ξ(h, ϕ(x 2 )), with ϕ(x 2 ) = γ 2 ϕ m f (x 2 ), (31)
where Ξ(h, ϕ) denotes a jump Markov process, occurring at a rate ϕ, whose amplitude is distributed with density h as given in (30). Analogously, in the Case 2 effector dominated situation the single cell equation becomes

dx 3 dt = -γ 3 x 3 + Ξ(h, ϕ(x 3 )), with ϕ(x 3 ) = γ 3 ϕ m f (x 3 ). ( 32
)
Equations 31 and 32 can both be written as

dx dt = -γx + Ξ(h, ϕ(x)), with ϕ(x) = γκ b f (x), κ b ≡ ϕ m .
Remark 4. In the case of bursting we will always take κ b ≡ ϕ m in contrast to the deterministic case where

κ d = b d ϕ m .
From Mackey and Tyran-Kamińska (2008) the corresponding operator equation for the evolution of the density u(t, x) when there is a single dominant slow variable is given by

∂u(t, x) ∂t -γ ∂(xu(t, x)) ∂x = -γκ b f (x)u(t, x) + γκ b x 0 f (y)u(t, y)h(x -y)dy. ( 33
)
Remark 5. This is a straightforward generalization of what Gardiner (1983, Section 3.4) refers to as the differential Chapman-Kolmogorov equation.

Stationary solutions u * (x) of ( 33) are solutions of

- d(xu * (x)) dx = -κ b f (x)u * (x) +κ b x 0 f (y)u * (y)h(x-y)dy. ( 34
)
If there is a unique stationary density, then the solution u(t, x) of Equation 33 is said to be asymptotically stable [START_REF] Lasota | Chaos, fractals, and noise[END_REF] in the sense that

lim t→∞ ∞ 0 |u(t, x) -u * (x)|dx = 0
for all initial densities u(0, x).

Theorem 3. (Mackey and Tyran-Kamińska, 2008, Theorem 7).

The unique stationary density of Equation 34, with f given by Equation 9and h given by ( 29), is

u * (x) = C x e -x/b exp κ b x f (y) y dy , ( 35
)
where C is a normalizing constant such that ∞ 0 u * (x)dx = 1. Further, u(t, x) is asymptotically stable. Remark 6. The stationary density ( 35) is found by rewriting Equation 34 in the form

dy(x) dx + y(x) b -κ b f (x) x y(x) = 0, y(x) ≡ xu * (x)
using Laplace transforms and solving by quadratures. Note also that we can represent u * as

u * (x) = C exp x κ b f (y) y - 1 b - 1 y dy,
where C is a normalizing constant.

Distributions in the presence of bursting 4.2.1. Protein distribution in the absence of control

If the burst frequency ϕ = γκ b f is independent of the level of all of the participating molecular species, then the solution given in Equation 35 is the density of the gamma distribution:

u * (x) = 1 b κ b Γ(κ b ) x κ b -1 e -x/b ,
where Γ(•) denotes the gamma function. For κ b ∈ (0, 1), u * (0) = ∞ and u * is decreasing while for κ b > 1, u * (0) = 0 and there is a maximum at x = b(κ b -1).

Controlled bursting

We next consider the situation in which the burst frequency ϕ is dependent on the level of x, c.f. Equation 5. This requires that we evaluate κ b

x f (y)

y dy = x κ b y 1 + y n Λ + ∆y n dy = ln x κ b Λ -1 (Λ + ∆x n ) θ ,
where Λ, ∆ are enumerated in Table 1 for both the inducible and repressible operons treated in Section 2.2 and

θ = κ b n∆ 1 - ∆ Λ .
Consequently, the steady state density (35) explicitly becomes The first two terms of Equation 36 are simply proportional to the density of the gamma distribution. For 0 < κ b Λ -1 < 1 we have u * (0) = ∞ while for κ b Λ -1 > 1, u * (0) = 0 and there is at least one maximum at a value of x > 0. We have u * (x) > 0 for all x > 0 and from Remark 6 it follows that

u * (x) = Ce -x/b x κ b Λ -1 -1 (Λ + ∆x n ) θ . ( 36
u ′ * (x) = u * (x) κ b f (x) x - 1 b - 1 x , x > 0. ( 37
)
Observe that if κ b ≤ 1 then u * is a monotone decreasing function of x, since κ b f (x) ≤ 1 for all x > 0. Thus we assume in what follows that κ b > 1.

Since the analysis of the qualitative nature of the stationary density leads to different conclusions for the inducible and repressible operon cases, we consider each in turn.

Bursting in the inducible operon

For θ > 0, as in the case of an inducible operon, the third term of Equation 36 is a monotone increasing function of x and, consequently, there is the possibility that u * may have more than one maximum, indicative of the existence of bistable behaviour. In this case, the stationary density becomes

u * (x) = Ce -x/b x κ b K -1 -1 (K + x n ) θ , θ = κ b n (1 -K -1 ).
From (37) it follows that we have u ′ * (x) = 0 for x > 0 if and only if 1

κ b x b + 1 = 1 + x n K + x n . ( 38
)
Again, graphical arguments (see Figure 5) show that there may be up to three roots of (38). For illustrative values of n, K, and b, Figure 6 shows the graph of the values of x at which u ′ * (x) = 0 as a function of κ b . When there are three roots of (38), we label them as x1 < x2 < x3 . Generally we cannot determine when there are three roots. However, we can determine when there are only two roots x1 < x3 from the argument of Section 2.3.2. At x1 and x3 we will not only have Equation 38 satisfied but the graph of the right hand side of (38) will be tangent to the graph of the left hand side at one of them so the slopes will be equal. Differentiation of (38) yields the second condition

n x n-1 (K + x n ) 2 = 1 κ b b(K -1) (39)
We first show that there is an open set of parameters (b, K, κ b ) for which the stationary density u * is bimodal. From Equations 38 and 39 it follows that the value of x ± at which tangency will occur is given by

x ± = b(κ b -1)z ±
and z ± are positive solutions of equation

z n = 1 -z -β(1 -z) 2 , where β = K(κ b -1) (K -1)κ b .
We explicitly have

z ± = 1 2βn 2βn -(n + 1) ± (n + 1) 2 -4βn provided that (n + 1) 2 4n ≥ β = K(κ b -1) (K -1)κ b . ( 40
)
Equation 40 is always satisfied when κ b < K or when κ b > K and K is as in the deterministic case (19). Observe also that we have z + > 0 > z -for κ b < K and z + > z -> 0 for κ b > K. The two corresponding values of b at which a tangency occurs are given by We now want to find the analogy between the bistable behavior in the deterministic system and the existence of bimodal stationary density u * . To this end we fix the parameters b > 0 and K > 1 and vary κ b as in Figure 5. Equations 38 and 39 can also be combined to give an implicit equation for the value of x ± at which tangency will occur

b ± = 1 (κ b -1)z ± n K β(1 -z ± ) -K and z ± > 0. If κ b < K then u * (0) = ∞
x 2n -(K -1) n - K + 1 K -1 x n -nb(K -1)x n-1 + K = 0 (41)
and the corresponding values of κ b± are given by

κ b± = x ∓ + b b K + x n ∓ 1 + x n ∓ . ( 42
)
There are two cases to distinguish. Case 1. 0 < κ b < K. In this case, u * (0) = ∞. Further, the same graphical considerations as in the deterministic case show that there can be none, one, or two positive solutions to Equation 38 Remark 7. Remember that the case n = 1 cannot display bistability in the deterministic case. However, in the case of bursting in the inducible system when n = 1, if K b + 1 < κ b < K and b > K K-1 , then u * (0) = ∞ and u * also has a maximum at x3 > 0. Thus in this case one can have a Bimodal type 1 stationary density.

We now choose to see how the average burst size b affects bistability in the density u * by looking at the parametric plot of

κ b (x) versus K(x). Define F(x, b) = x n + 1 nx n-1 (x + b) . (43)
Then

K(x, b) = 1 + x n F(x, b) 1 -F(x, b) and κ b (x, b) = [K(x, b)+x n ] x + b b(x n + 1)
.

(44) The bifurcation diagram obtained from a parametric plot of K versus κ b (with x as the parameter) is illustrated in Figure 7 for n = 4 and two values of b. Note that it is necessary for 0 < K < κ b in order to obtain Bimodal type 2 behaviour.

For bursting behaviour in an inducible situation, there are two different bifurcation patterns that are possible. The two different cases are delineated by the respective values of K and κ b , as shown in Figure 6 and Figure 7. Both bifurcation scenarios share the property that while increasing the bifurcation parameter κ b from 0 to ∞, the stationary density u * passes from a 10 . Though somewhat obscured by the logarithmic scale for x, the graphs always intersect the κ b axis at κ b = K. Additionally, it is important to note that u ′ * (0) = 0 for K < κ b , and that there is always a maximum at 0 for 0 < κ b < K. See the text for further details. unimodal density with a peak at a low value (either 0 or x1 ) to a bimodal density and then back to a unimodal density with a peak at a high value ( x3 ). In what will be referred as Bifurcation type 1, the maximum at x = 0 disappears when there is a second peak at x = x3 . The sequence of densities encountered for increasing values of κ b is then: Unimodal type 1 to a Bimodal type 1 to a Bimodal type 2 and finally to a Unimodal type 2 density.

In the Bifurcation type 2 situation, the sequence of density types for increasing values of κ b is: Unimodal type 1 to a Unimodal type 2 and then a Bimodal type 2 ending in a Unimodal type 2 density.

The two different kinds of bifurcation that can occur are easily illustrated for b = 1 as the parameter κ b is increased. (An enlarged diagram in the region of interest is shown in Figure 8.) Figure 9 illustrates Bifurcation type 1, when K = 4, and κ b increases from low to high values. As κ b increases, we pass from a Unimodal type 1 density, to a Bimodal type 1 density. Further increases in κ b lead to a Bimodal type 2 density and finally to a Unimodal type 2 density. This bifurcation cannot occur, for example, when b = 1 10 and K ≤ 15 (see Figure 7). Figure 10 shows Bifurcation type 2, when K = 3. As κ b increases, we pass from a Unimodal type 1 density, to a Unimodal type 2 density. Then with further increases in κ b , we pass to a Bimodal type 2 density and finally back to a Unimodal type 2 density.

Remark 8. There are several qualitative conclusions to be drawn from the analysis of this section.

1. The presence of bursting can drastically alter the regions of parameter space in which bistability can occur relative to the deterministic case. 5.9, u * has two local maxima at 0 < x1 < x3 . Finally, for κ b 5.9, u * has a single maximum at x3 > 1. Note that for each plot of the density, the scale of the ordinate is arbitrary to improve the visualization.

2. When 0 < κ b < K, at a fixed value of κ b , increasing the average burst size b can lead to a bifurcation from Unimodal type 1 to Bimodal type 1. 3. When 0 < K < κ b , at a fixed value of κ b , increasing b can lead to a bifurcation from Unimodal type 2 to Bimodal type 2 and then back to Unimodal type 2.

Bursting in the repressible operon

The possible behaviours in the stationary density u * for the repressible operon are easy to delineate based on the analysis of the previous section, with Equation 38 replaced by

1 κ b x b + 1 = 1 + x n 1 + ∆x n . (45)
Again graphical arguments (see Figure 12) show that Equation 45 may have either none or one solution. Namely, 1. For 0 < κ b < 1, u * (0) = ∞ and u * is decreasing. Equation 45does not have any solution (Unimodal type 1). 2. For 1 < κ b , u * (0) = 0 and u * has a single maximum at a value of x > 0 determined by the single positive solution of Equation 45 (Unimodal type 2).

Recovering the deterministic case

We can recover the deterministic behaviour from the bursting dynamics with a suitable scaling of the parameters and limiting procedure. With bursting production there are two important parameters (the frequency κ b and the amplitude b), while with deterministic production there is only κ d . The natural limit to consider is when 4.45 u * has a single maximum at x3 > 0. Note that for each plot of the density, the scale of the ordinate is abritrary to improve the visualization.

b → 0, κ b → ∞ with bκ b ≡ κ d .
In this limit, the implicit equations which define the maximum points of the steady state density, become the implicit equations ( 15) and ( 16) which define the stable steady states in the deterministic case.

The bifurcations will also take place at the same points, because we recover Equation 18 in the limit. However, Bimodality type 1 as well as the Unimodal type 1 behaviours will no longer be present, as in the deterministic case, because for κ b → ∞ we have κ b > K. Finally, from the analytical expression for the steady-state density (36) u * will became more sharply peaked as b → 0. Due to the normalization constant (which depends on b and κ b ), the mass will be more concentrated around the larger maximum of u * .

Distributions with fluctuations in the degradation rate

Generalities

For a generic one dimensional stochastic differential equation of the form

dx(t) = α(x)dt + σ(x)dw(t)
the corresponding Fokker Planck equation the problem. In the particular case when J = 0 at one of the boundaries (a reflecting boundary) then J = 0 for all x in the domain and the steady state solution u * of Equation 46 is easily obtained with a single quadrature as

∂u ∂t = - ∂(αu) ∂x + 1 2 ∂ 2 (σ 2 u) ∂x 2 (
u * (x) = C σ 2 (x) exp 2 x α(y) σ 2 (y) dy ,
where C is a normalizing constant as before.

Fluctuations in degradation rate

In our considerations of the effects of continuous fluctuations, we examine the situation in which fluctuations appear in the degradation rate γ of the generic equation (28). From standard chemical kinetic arguments [START_REF] Oppenheim | Stochastic and deterministic formulation of chemical rate equations[END_REF], if the fluctuations are Gaussian distributed the mean numbers of molecules decaying in a time dt is simply γxdt and the standard deviation of these numbers is proportional to √ x. Thus we take the decay to be given by the sum of a deterministic component γxdt and a stochastic component σ √ xdw(t), where w is a standard Brownian motion, and write Equation 28 as a stochastic differential equation in the form

dx = γ[κ d f (x) -x]dt + σ √ xdw.
Within the Ito interpretation of stochastic integration, this equation has a corresponding Fokker Planck equation for the evolution of the ensemble density u(t, x) given by [START_REF] Lasota | Chaos, fractals, and noise[END_REF]) In the situation we consider here, σ(x) = σ √ x and α(x) = γκ d f (x) -γx. Further, since concentrations of molecules cannot become negative the boundary at x = 0 is reflecting and the stationary solution of Equation 47 is given by

∂u ∂t = - ∂ (γκ d f (x) -γx)u ∂x + σ 2 2 ∂ 2 (xu) ∂x 2 . ( 47
u * (x) = C x e -2γx/σ 2 exp 2γκ d σ 2
x f (y) y dy .

Set κ e = 2γκ d /σ 2 . Then the steady state solution is given explicitly by

u * (x) = Ce -2γx/σ 2 x κ e Λ -1 -1 [Λ + ∆x n ] θ , (48)
where Λ, ∆ ≥ 0 and θ are given in Table 1.

Remark 9. Two comments are in order.

1. Because the form of the solutions for the situation with bursting (intrinsic noise) and extrinsic noise are identical, all of the results of the previous section can be carried over here with the proviso that one replaces the average burst amplitude b with b → σ 2 /2γ ≡ b w and κ b → κ e = 2γκ d /σ 2 ≡ κ d /b w . 2. We can look for the regions of bimodality in the (K, κ d )plane, for a fixed value of b w . We have the implicit equation for x ±

x 2n -(K -1) n - K + 1 K -1 x n -nb w (K -1)x n-1 + K = 0
and the corresponding values of κ d are given by

κ d± = (x ∓ + b w ) K + x n ∓ 1 + x n ∓ .
Then the bimodality region in the (K, κ d )-plane with noise in the degradation rate is the same as the bimodality region for bursting in the (K, bκ b )-plane.

We have also the following result.

Theorem 4. (Pichór and Rudnicki, 2000, Theorem 2). The unique stationary density of Equation 47 is given by Equation 48. Further u(t, x) is asymptotically stable.

The deterministic limit

Here again we can recover the deterministic behavior from a limit in the extrinsic fluctuations dynamics. In this case, however, the frequency and the amplitude of the perturbation are already scaled. Then the limit σ → 0 gives the same result as in the deterministic case.

Discussion and conclusions

In trying to understand experimentally observed distributions of intracellular components from a modeling perspective, the norm in computational and systems biology is often to use algorithms developed initially by [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF] to solve the chemical master equation for specific situations. See [START_REF] Lipniacki | Transcriptional stochasticity in gene expression[END_REF] for a typical example. However these investigations demand long computer runs, are computationally expensive, and further offer little insight into the possible diversity of behaviours that different gene regulatory networks are capable of.

There have been notable exceptions in which the problem has been treated from an analytical point of view, c.f. [START_REF] Kepler | Stochasticity in transcriptional regulation: Origins, consequences, and mathematical representations[END_REF], [START_REF] Friedman | Linking stochastic dynamics to population distribution: An analytical framework of gene expression[END_REF], [START_REF] Bobrowski | Asymptotic behavior of distributions of mRNA and protein levels in a model of stochastic gene expression[END_REF], and Shahrezaei and Swain (2008a). The advantage of an analytic development is that one can determine how different elements of the dynamics shape temporal and steady state results for the densities u(t, x) and u * (x) respectively.

Here we have extended this analytic treatment to simple situations in which there is either bursting transcription and/or translation (building on and expanding the original work of [START_REF] Friedman | Linking stochastic dynamics to population distribution: An analytical framework of gene expression[END_REF]), or fluctuations in degradation rates, as an alternative to the [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF] algorithm approach. The advantage of the analytic approach that we have taken is that it is possible, in some circumstances, to give precise conditions on the statistical stability of various dynamics. Even when analytic solutions are not available for the partial integro-differential equations governing the density evolution, the numerical solution of these equations may be computationally more tractable than using the [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF] approach.

One of the more surprising results of the work reported here is that the stationary densities in the presence of bursting noise derived in Section 4 are analytically indistinguishable from those in the presence of degradation noise studied in Section 5. We had expected that there would be clear differences that would offer some guidance for the interpretation of experimental data to determine whether one or the other was of predominant importance. Of course, the next obvious step is to examine the problem in the presence of both noise sources simultaneously. However, the derivation of the evolution equation in this case, as has been pointed out [START_REF] Hierro | Singular boundaries in the forward Chapman-Kolmogorov differential equation[END_REF], is not straightforward and we will report on our results in a separate communication.

In terms of the issue of when bistability, or a unimodal versus bimodal stationary density is to be expected, we have pointed out the analogy between the unimodal and bistable behaviours in the deterministic system and the existence of bimodal stationary densities in the stochastic systems. Our analysis makes clear the critical role of the dimensionless parameters n, κ (be it κ d , κ b , or κ e ), b (either b or b w ), and the fractional leakage K -1 . The relations between these defining the various possible behaviours are subtle, and we have given these in the relevant sections of our analysis.

The appearance of both unimodal and bimodal distributions of molecular constituents as well as what we have termed Bifurcation Type 1 and Bifurcation Type 2 have been extensively discussed in the applied mathematics literature (c.f. [START_REF] Horsthemke | Noise Induced Transitions: Theory and Applications in Physics, Chemistry, and Biology[END_REF], [START_REF] Feistel | Evolution of Complex Systems[END_REF] and others) and the bare foundations of a stochastic bifurcation theory have been laid down by [START_REF] Arnold | Random dynamical systems[END_REF]. Significantly, these are also well documented in the experimental literature as has been shown by [START_REF] Gardner | Construction of a genetic toggle switch in Escherichia coli[END_REF], [START_REF] Acar | Enhancement of cellular memory by reducing stochastic transitions[END_REF], [START_REF] Friedman | Linking stochastic dynamics to population distribution: An analytical framework of gene expression[END_REF], [START_REF] Hawkins | The regulatory roles of the galactose permease and kinase in the induction response of the GAL network in Saccharomyces cerevisiae[END_REF], [START_REF] Zacharioudakis | A yeast catabolic enzyme controls transcriptional memory[END_REF], [START_REF] Mariani | Short-term memory in gene induction reveals the regulatory principle behind stochastic IL-4 expression[END_REF][START_REF] Song | Estimating the stochastic bifurcation structure of cellular networks[END_REF] for both prokaryotes and eukaryotes. If the biochemical details of a particular system are sufficiently well characterized from a quantitative point of view so that relevant parameters can be estimated, it may be possible to discriminate between whether these behaviours are due to the presence of bursting transcription/translation or extrinsic noise.

Figure 1 :

 1 Figure 1: Schematic illustration of the possibility of one, two or three solutions of Equation 15 for varying values of κ d with inducible regulation. The monotone increasing graph is the function f of Equation 10, and the straight lines correspond to x/κ d for (in a clockwise direction) κ d ∈ [0, κ d-), κ d = κ d-, κ d ∈ (κ d-, κ d+ ), κ d = κ d+ , and κ d+ < κ d . This figure was constructed with n = 4 and K = 10 for which κ d-= 3.01 and κ d+ = 5.91 as computed from (18). See the text for further details.

Figure 2 :

 2 Figure 2: Full logarithmic plot of the steady state values of x * versus κ d for an inducible system, obtained from Equation 15, for n = 4 and K = 2, 5, 10, and 15 (left to right) illustrating the dependence of the occurrence of bistability on K. See the text for details.

Figure 3 :

 3 Figure3: In this figure we present a parametric plot (for n = 4) of the bifurcation diagram in (K, κ d ) parameter space delineating one from three steady states in a deterministic inducible operon as obtained from Equations 15 and 16. The upper (lower) branch corresponds to κ d-(κ d+ ), and for all values of (K, κ d ) in the interior of the cone there are two locally stable steady states X * 1 , X * 3 , while outside there is only one. The tip of the cone occurs at (K, κ d ) = ((5/3) 2 , (5/3) 4 √ 5/3) as given by Equations 19 and 20. For K ∈ [0, (5/3) 2 ) there is but a single steady state.

Figure 4 :

 4 Figure 4: Schematic illustration that there is only a single solution of Equation 15 for all values of κ d with repressible regulation. The monotone decreasing graph is f for a repressible operon, while the straight lines are x/κ d . This figure was constructed with n = 4 and ∆ = 10. See the text for further details.

Figure 5 :

 5 Figure 5: Schematic illustration of the possibility of one, two or three solutions of Equation 38 for varying values of κ b with bursting inducible regulation. The straight lines correspond (in a clockwise direction) to κ b ∈ (0, κ b-), κ b = κ b-, κ b ∈ (κ b-, κ b+ ) (and respectively κ b < K, κ b = K, K < κ b ), κ b = κ b+ , and κ b+ < κ b . This figure was constructed with n = 4, K = 10 and b = 1 for which κ b-= 4.29 and κ b+ = 14.35 as computed from (42). See the text for further details.

  and u * is decreasing for b ≤ b + , while for b > b + there is a local maximum at x > 0. If κ b > K then u * (0) = 0 and u * has one or two local maximum. As a consequence, for n > 1 we have a bimodal steady state density u * if and only if the parameters κ b and K satisfy (40), κ b > K, and b ∈ (b + , b -).

  . If κ b < κ b-, there are no positive solutions, u * is a monotone decreasing function of x. If κ b > κ b-, there are two positive solutions ( x2 and x3 in our previous notation, x1 has become negative and not of importance) and there will be a maximum in u * at x3 with a minimum in u * at x2 . Case 2. 0 < K < κ b . Now, u * (0) = 0 and there may be one, two, or three positive roots of Equation38. We are interested in knowing when there are three which we label as x1 < x2 < x3 as x1 , x3 will correspond to the location of maxima in u * while x2 will be the location of the minimum between them and the condition for the existence of three roots is κ b-< κ b < κ b+ .We see then that the different possibilities depend on the respective values of K, κ b-, κ b+ , and κ b . To summarize, we may characterize the stationary density u * for an inducible operon in the following way:1. Unimodal type 1: u * (0) = ∞ and u * is decreasing for 0 < κ b < κ b-and 0 < κ b < K 2.Unimodal type 2: u * (0) = 0 and u * has a single maximum at (a) x1 > 0 for K < κ b < κ b-or (b) at x3 > 0 for κ b+ < κ b and K < κ b 3. Bimodal type 1: u * (0) = ∞ and u * has a single maximum at x3 > 0 for κ b-< κ b < K 4. Bimodal type 2: u * (0) = 0 and u * has two maxima at x1 , x3 , 0 < x1 < x3 for κ b-< κ b < κ b+ and K < κ b

Figure 6 :

 6 Figure 6: Full logarithmic plot of the values of x at which u ′ * (x) = 0 versus the parameter κ b , obtained from Equation 38, for n = 4, K = 10, and (left to right) b = 5, 1 and b = 110 . Though somewhat obscured by the logarithmic scale for x, the graphs always intersect the κ b axis at κ b = K. Additionally, it is important to note that u ′ * (0) = 0 for K < κ b , and that there is always a maximum at 0 for 0 < κ b < K. See the text for further details.

Figure 7 :Figure 8 :

 78 Figure 7: In this figure we present two bifurcation diagrams (for n = 4) in (K, κ b ) parameter space delineating unimodal from bimodal stationary densities u * in an inducible operon with bursting as obtained from Equations 44 with 43. The upper cone-shaped plot is for b = 1 10 while the bottom one is for b = 1. In both cone shaped regions, for any situation in which the lower branch is above the line κ b = K (lower straight line) then bimodal behaviour in the stationary solution u * (x) will be observed with maxima in u * at positive values of x, x1 and x3 .

Figure 9 :

 9 Figure 9: In this figure we illustrate Bifurcation type 1 when intrinsic bursting is present. For a variety of values of the bifurcation parameter κ b (between 3 and 6 from top to down), the stationary density u * is plotted versus x between 0 and 8. The values of the parameters used in this figure are b = 1, K = 4, and n = 4. For κ b 3.5, u * has a single maximum at x = 0. For 3.5 κ b < 4, u * has two local maxima at x = 0 and x3 > 1. For 4 < κ b5.9, u * has two local maxima at 0 < x1 < x3 . Finally, for κ b 5.9, u * has a single maximum at x3 > 1. Note that for each plot of the density, the scale of the ordinate is arbitrary to improve the visualization.

Figure 10 :

 10 Figure 10: An illustration of Bifurcation type 2 for intrinsic bursting. For several values of the bifurcation parameter κ b (between 2.8 and 5 from top to down), the stationary density u * is plotted versus x between 0 and 8. The parameters used are b = 1, K = 3, and n = 4. For κ b < 3, u * has a single maximum at x = 0, and for 3 < κ b 3.3, u * has a single maximum at x1 > 0. For 3.3 κ b 4.45, u * has two local maxima at 0 < x1 < x3 , and finally for κ b4.45 u * has a single maximum at x3 > 0. Note that for each plot of the density, the scale of the ordinate is abritrary to improve the visualization.

Figure 11 :

 11 Figure 11: The presence of bursting can drastically alter regions of bimodal behaviour as shown in this parametric plot (for n = 4) of the boundary in (K, b • κ b ) parameter space delineating unimodal from bimodal stationary densities u * in an inducible operon with bursting and in (K, κ d ) parameter space delineating one from three steady states in the deterministic inducible operon. From top to bottom, the regions are for b = 10, b = 1, b = 0.1 and b = 0.01. The lowest (heavy dashed line) is for the deterministic case. Note that for b = 0.1, the two regions of bistability and bimodality coincide and are indistinguishable from one another.

Figure 12 :

 12 Figure 12: Schematic illustration that there can be one or no solution of Equation 45, depending on the value of κ b , with repressible regulation. The straight lines correspond (in a clockwise direction) to κ b = 2 and κ b = 0.8. This figure was constructed with n = 4, ∆ = 10 and b = 1. See the text for further details.

Table 1 :

 1 Definitions of the parameters A, B, Λ, ∆ and θ. See the text and Section 2.2 for more detail.

Acknowledgements

This work was supported by the Natural Sciences and Engineering Research Council (NSERC, Canada), the Mathematics of Information Technology and Complex Systems (MITACS, Canada), the Alexander von Humboldt Stiftung, the State Committee for Scientific Research (Poland) and the Ecole Normale Superieure Lyon (ENS Lyon, France). The work was carried out at McGill University, Silesian University, the University of Bremen, and the Oxford Centre for Industrial and Applied Mathematics (OCIAM), University of Oxford.


Agostino Martinelli
email: agostino.martinelli@inria.fr

Overdamped 2D-Brownian motion for self-propelled and nonholonomic particles

This paper provides a complete analytical solution to the Brownian motion of a self-propelled particle that moves in two spatial dimensions in the overdamped limit by satisfying kinematic constraints. From a mathematical point of view, it provides the analytical solution of a Fokker Plank equation for which the detailed balance is not verified. The analytical expression of any-order moment of the probability distribution is obtained by a direct integration of the Langevin equation. Each moment is expressed as a multiple integral of the active motion performed by the particle. For special active motions, these integrals can be easily solved. In particular, they can be expressed as the real or the imaginary part of suitable analytic functions. This is explicitly shown for the first two moments.

Introduction

In recent years a great interest has been devoted to investigate the statistical properties of the motion of active (or self-propelled) particles. These particles differ from passive particles since they move under the action of an internal force. Typical examples of these particles are Brownian motors [START_REF] Mark | On Diffusion in Narrow Random Channels[END_REF][START_REF] Peter | Artificial Brownian motors: Controlling transport on the nanoscale[END_REF][START_REF] Peter | Brownian motors: noisy transport far from equilibrium[END_REF], biological and artificial microswimmers [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF][START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF][START_REF] Leoni | A basic swimmer at low Reynolds number[END_REF][START_REF] Sokolov | Concentration Dependence of the Collective Dynamics of Swimming Bacteria[END_REF][START_REF] Tierno | Magnetically actuated colloidal microswimmers[END_REF][START_REF] Tierno | Controlled swimming in confined fluids of magnetically actuated colloidal rotors[END_REF][START_REF] Trouilloud | Soft swimming: Exploiting deformable interfaces for low-Reynolds number locomotion[END_REF], macroscopic animals [START_REF] Bazazi | Nutritional state and collective motion: from individuals to mass migration[END_REF][START_REF] Couzin | Collective memory and spatial sorting in animal groups[END_REF][START_REF] Niwa | Self-organizing dynamic model offish schooling[END_REF] and even pedestrians [START_REF] Turner | Brownian Motion Applied to Human Intersections[END_REF]. As in the case of passive particles, the motion can often be characterized by the Langevin equation in the overdamped limit [START_REF] Van Teeffelen | Dynamics of a Brownian circle swimmer[END_REF]. In this case, the probability distribution satisfies the Smoluchowski equation, which is a special case of the Fokker-Plank equation [START_REF] Risken | The Fokker-Planck Equation Methods of Solution and Applications[END_REF]. On the other hand, since in most of cases of active particles the so-called detailed balance condition with respect to a given time reversal transformation is not satisfied, it is more convenient to derive the properties of the probability distribution by using the stochastic differential equation (i.e., the Langevin equation) instead of directly trying to solve the Fokker-Plank equation (we remind the reader that analytic solutions of the Fokker-Plank equation can be found when the detailed balance condition is satisfied: when this condition is not satisfied, it is in general very difficult to find analytic solutions in the non linear case and for a large number of variables [START_REF] Risken | The Fokker-Planck Equation Methods of Solution and Applications[END_REF]). Very recently, the Langevin equation has successfully been used to analytically compute the moments of the probability distribution for active particles [6,7,[START_REF] Van Teeffelen | Dynamics of a Brownian circle swimmer[END_REF]. Specifically, the first two moments have been derived for rodlike particles [START_REF] Van Teeffelen | Dynamics of a Brownian circle swimmer[END_REF] and spherical particles [6]. In both these cases the orientation of the particle was confined to two dimensions. The work presented in [7] extends these works by including the three-dimensional case for both isotropic and anisotropic particles. Additionally, the first four moments were analytically derived (instead of the first two).

In this paper we consider the case of active particles that satisfy nonholonomic constraints. Specifically, our stochastic motion model (presented in section 2) differs from the one in [7] since both the active motion and the noise must satisfy the kinematic constraint. Additionally, the active motion is time-dependent and it is characterized by two independent functions of time, the former characterizes the linear speed, the latter the torque. On the other hand, we restrict the computation to the two-dimensional case even if we believe that the extension to the three-dimensional case is possible. Starting from this motion model, we compute the analytical expression of any-order moment of the probability distribution (section 3). Specifically, each moment is expressed as a multiple integral of the active motion performed by the particle. In other words, the expression is a multiple integral on the two time-dependent functions previously mentioned. Finally, when the ratio of these two functions is independent of time, these integrals can be easily solved. In particular, the active trajectory and the noise can be characterized by a suitable complex quantity and the integrals can be expressed as the real or the imaginary part of suitable analytic functions of this complex quantity. This is explicitly shown for the first two moments in section 4. Conclusions are provided in section 5.

Note that we already computed the statistics up to the second order for a similar motion model in the framework of mobile robotics [START_REF] Martinelli | The odometry error of a mobile robot with a synchronous drive system[END_REF][START_REF] Martinelli | The accuracy on the parameter estimation of an odometry system of a mobile robot[END_REF]. (in the expressions in [START_REF] Martinelli | The odometry error of a mobile robot with a synchronous drive system[END_REF] there are some typos, corrected in [START_REF] Martinelli | The accuracy on the parameter estimation of an odometry system of a mobile robot[END_REF]). Here, we extend the computation in [START_REF] Martinelli | The odometry error of a mobile robot with a synchronous drive system[END_REF] by including any-order statistics. As it will be seen, this extension requires considerable analytical effort. We also show that the expressions hold for a much more general model (the one described in section 2) and not only for a simple specific odometry error model (as shown in [START_REF] Martinelli | The odometry error of a mobile robot with a synchronous drive system[END_REF]). As a result, this paper provides a complete analytical solution to a Fokker-Plank equation for which the detailed balance condition is not verified ‡.

Stochastic motion model

We consider a particle that moves in a 2D-environment. The configuration of the particle is characterized by its position and orientation. In Cartesian coordinates, we denote them by r ≡ [x, y] T and θ, respectively. We assume that the motion of the particle satisfies the unicycle model [START_REF] Latombe | Robot motion planning[END_REF], which is the most general nonholonomic model for 2D-motions, where the shift of the particle can only occur along one direction ([cos θ, sin θ] T ):

    ẋ = v cos θ ẏ = v sin θ θ = ω (1)
The quantities v = v(t) and ω = ω(t) are two functions of time and are the linear and angular speed, respectively. We assume that these functions are known and we call the motion that results from them the active motion. Now we want to introduce a stochastic model that generalizes equation (1) by accounting a Gaussian white noise. We assume that also the noise satisfies the same nonholonomic constraint. The stochastic differential equations are:

    dx(t) = cos θ(t) [v(t)dt + α(t)dw r (t)] dy(t) = sin θ(t) [v(t)dt + α(t)dw r (t)] dθ(t) = ω(t)dt + β(t)dw θ (t) (2)
where [dw r (t), dw θ (t)] T is a standard Wiener process of dimension two [START_REF] Oksendal | Stochastic Differential Equations: An Introduction with Applications[END_REF]. The functions α(t) and β(t) are modelled according to the following physical requirement (diffusion in the overdamped regime). We require that, when the particle moves during the infinitesimal time interval dt, the shift is a random Gaussian variable, whose variance increases linearly with the travelled distance. This is obtained by setting α(t) = K r |v(t)|, with K r a positive parameter that characterizes our system (interaction particle-environment). Similarly, we require that the rotation accomplished during the same time interval dt is a random Gaussian variable, whose variance increases linearly with the travelled distance. Hence, we set β(t) = K θ |v(t)|, with K θ another positive parameter that characterizes our system (interaction particle-environment). The equations in (2) are the Langevin equations in the overdamped limit. From now on, we assume that the particle can only move ahead, i.e., v(t) ≥ 0. Instead of the time t, we use the curve length ds = v(t)dt. Note that v(t) is a deterministic and known function of time. Additionally, we will denote by µ the ratio between the angular and the linear speed. Equations in (2) read:

     dx(s) = cos θ(s) [ds + K 1 2 r dw r (s)] dy(s) = sin θ(s) [ds + K 1 2 r dw r (s)] dθ(s) = µ(s)ds + K 1 2 θ dw θ (s) (3)
The associated Smoluchowski equation is:

∂p ∂s = -∇ • (D 1 p) + ∇ • (D 2 ∇p) = 0 (4)
where:

• p = p(x, y, θ; s(t)) is the probability density for the particle at the configuration (x, y, θ) and at time t;

• D 1 = [cos θ, sin θ, µ] T is the drift vector;

• D 2 = 1 2     K r cos 2 θ K r cos θ sin θ 0 K r cos θ sin θ K r sin 2 θ 0 0 0 K θ   
 is the diffusion tensor.

Our goal is to obtain the probability density p(x, y, θ; s). This will be done in the next two sections.

Computation of the probability distribution

The probability density p(x, y, θ; s) satisfies the Smoluchowski partial differential equation in ( 4), which is a special case of the Fokker-Plank equation. Since the detailed balance condition is not satisfied [START_REF] Risken | The Fokker-Planck Equation Methods of Solution and Applications[END_REF], we follow a different procedure to compute p(x, y, θ; s). Specifically, we use the Langevin equation in (3) to compute the moments, up to any order, of the probability distribution. First of all, we remark that the third equation in ( 3) is independent of the first two, is independent of θ and is linear in dw θ . As a result, the probability density only in terms of θ, i.e., the probability P θ (θ; s) ≡ dx dy p(x, y, θ; s), is a Gaussian distribution with mean value θ 0 + s 0 ds ′ µ(s ′ ) and variance K θ s, where θ 0 is the initial orientation. This same result could also be obtained by integrating (4) in x and y and by using the divergence theorem. In other words, we have the following distribution for the orientation at a given value of s:

θ(s) = N (θ(s), K θ s) (5)
where N (•, •) denotes the normal distribution with mean value and variance the first and the second argument; θ(s) ≡ θ 0 + s 0 ds ′ µ(s ′ ). Note that, according to the unicycle model, a trajectory is completely characterized by its starting point and by the orientation vs the curve length, i.e., by the function θ(s). In the following, we will call the deterministic function θ(s), the active trajectory.

First and second-order statistics

Let us consider the first equation in [START_REF] Couzin | Collective memory and spatial sorting in animal groups[END_REF]. We compute the expression of x(s) by a direct integration. We divide the interval (0, s) in N equal segments, δs ≡ s N . We have:

x(s) = lim N →∞ N j=1 (δs + ǫ j ) cos θ j (6)
where ǫ j is a random Gaussian variable satisfying ǫ j = 0, ǫ j ǫ k = δ jk K r δs for j, k = 1, • • • , N (δ jk is the Kronecker delta) and θ j ≡ θ(jδs). On the other hand, we have:

θ j = θ(jδs) + j m=1 δθ m ≡ θ j + ∆θ j (7)
where, according to our stochastic model in (3), δθ j is a random Gaussian variable satisfying δθ j = 0, δθ j ǫ k = 0 and δθ j δθ k = δ jk K θ δs for j, k = 1, • • • , N . As a result, ∆θ j ≡ j m=1 δθ m is also a random Gaussian variable satisfying ∆θ j = 0,

∆θ j ∆θ k = K θ jδs if j ≤ k K θ kδs if j > k .
Starting from ( 6) and (7) and by remarking that cos ∆θ j = e -k θ jδs 2 and sin ∆θ j = 0, it is easy to obtain the mean value of x(s). We have:

x(s) = lim N →∞ N j=1 (δs + ǫ j ) cos θ j = lim N →∞ N j=1 δs cos θ j e -k θ jδs 2 = (8) = s 0 ds ′ cos θ(s ′ ) e -K θ s ′ 2
Similarly, it is possible to obtain the mean value of y(s):

y(s) = s 0 ds ′ sin θ(s ′ ) e -K θ s ′ 2 (9)
Finally, in a similar manner, but with some more computation, we obtain all the secondorder moments:

x(s) 2 = s 0 ds ′ s-s ′ 0 ds ′′ e -K θ s ′′ 2 (1 + χ c (s ′ )) cos[θ(s ′ + s ′′ ) -θ(s ′ )]+ (10) -χ s (s ′ ) sin[θ(s ′ + s ′′ ) -θ(s ′ )] + K r 2 s + s 0 ds ′ χ c (s ′ ) y(s) 2 = s 0 ds ′ s-s ′ 0 ds ′′ e -K θ s ′′ 2 (1 -χ c (s ′ )) cos[θ(s ′ + s ′′ ) -θ(s ′ )]+ (11) +χ s (s ′ ) sin[θ(s ′ + s ′′ ) -θ(s ′ )] + K r 2 s - s 0 ds ′ χ c (s ′ ) x(s) y(s) = s 0 ds ′ s-s ′ 0 ds ′′ e -K θ s ′′ 2 χ s (s ′ ) cos[θ(s ′ + s ′′ ) -θ(s ′ )]+ (12) +χ c (s ′ ) sin[θ(s ′ + s ′′ ) -θ(s ′ )] + K r 2 s 0 ds ′ χ s (s ′ ) σ xθ (s) ≡ x(s) θ(s) -x(s) θ(s) = 2K θ ∂ y(s) ∂K θ ( 13
)
σ yθ (s) ≡ y(s) θ(s) -y(s) θ(s) = -2K θ ∂ x(s) ∂K θ ( 14
)
where

χ c (s ′ ) ≡ cos[2θ(s ′ )] e -2K θ s ′ and χ s (s ′ ) ≡ sin[2θ(s ′ )] e -2K θ s ′ .
Obtaining the expression of higher-order moments will demand more tricky computation and will be dealt separately, in the next subsection. Here we conclude by considering the quantity D(s) 2 ≡ x(s) 2 + y(s) 2 , which provides the time-evolution of the square of the distance of the particle from its initial position. From ( 10) and ( 11) we obtain a simple expression for its mean value:

D(s) 2 = K r s + 2 s 0 ds ′ s-s ′ 0 ds ′′ e -K θ s ′′ 2 cos[θ(s ′ + s ′′ ) -θ(s ′ )] (15)

Computation of any-order moment

We introduce the following two complex random quantities:

u(s) ≡ lim N →∞ N j=1 (δs + ǫ j )e iθ j ; v(s) ≡ lim N →∞ N j=1 (δs + ǫ j )e -iθ j (16)
From ( 6) and ( 16) it is immediate to realize that x(s) = u(s)+v(s)

2

. Similarly we have y(s) = u(s)-v(s)

2i

. Hence, in order to compute any-order moment that involves x(s) and y(s) it suffices to compute u(s) p v(s) q for any p, q ∈ N . The computation of this quantity requires several tricky steps, which are provided in appendix Appendix A. The key is to separate all the independent random quantities in order to compute their mean values. This is obtained by arranging all the sums in a suitable manner (see appendix Appendix A for all the details). We have:

u(s) p v(s) q = ⌊ p+q 2 ⌋ n=0 K n r 2n l=0 p 2n -l q l m 2n -l m l m m! ( 17
) (2n -l -m -1)!!(l -m -1)!! s m e i(p-q)θ 0 2n -l -m 2 ! l -m 2 !(p -2n + l)!(q -l)! c s 0 ds 1 s s 1 ds 2 • • • s s β-1 ds β exp    β-1 b=0 i(p -q + Φ b )[θ(s b+1 ) -θ(s b )] - (p -q + Φ b ) 2 (s b+1 -s b )K θ 2   
where:

• ⌊ p+q 2 ⌋ is the largest integer not greater than p+q 2 ; • the second sum on l (i.e., 2n l=0 ) is restricted to the values of l for which 2nl ≤ p and l ≤ q (or, equivalently, we are using the convention that x y = 0 when y > x); • the sum on m (i.e., m ) goes from 0 to the minimum between l and 2nl and it is restricted to the integers m with the same parity of l (hence, both 2n-l-m 2 and l-m 2 are integers); • the symbols "!" and "!!" denote the factorial and the double factorial, respectively [START_REF] Callan | A combinatorial survey of identities for the double factorial[END_REF] (note that 0! = 0!! = (-1)!! = 1);

• β = p + qnm is the dimension of the remaining multiple integral (note that a multiple integral of dimension m has already been computed and provided the term s m );

• when β = 0 the sum on c, i.e. c , must be replaced with 1;

• when β ≥ 1 the sum on c is the sum over all the vectors c of dimension β, whose entries are -2, -1, 1 and 2: specifically, each vector c has l-m 2 entries equal to 2, ql equal to 1, p -2n + l equal to -1 and 2n-l-m 2 equal to -2 (note that the sum

c consists of β q-l β-q+l p-2n+l n-m l-m 2
addends);

• Φ b ≡ b a=1 c a (note that Φ β = qp); • s 0 ≡ 0 and θ(s 0 ) = θ 0 (note that s 0 is not a variable of integration).

In order to complete the derivation of the statistics for our problem, we need to compute any-order moment that also involves the orientation θ. We provide the formula for the quantity u(s) p v(s) q θ(s) r , where θ(s) ≡ θ(s)θ(s). The details of this computation are provided in appendix Appendix B. We have, for any p, q, r ∈ N :

u(s) p v(s) q θ(s) r = ⌊ p+q 2 ⌋ n=0 K n r 2n l=0 p 2n -l q l m 2n -l m l m m! (18) (2n -l -m -1)!!(l -m -1)!! s m e i(p-q)θ 0 2n -l -m 2 ! l -m 2 ! (p -2n + l)!(q -l)! c γ r!(γ β+1 -1)!!K r 2 θ β b=0 γ b+1 ! s 0 ds 1 s s 1 ds 2 • • • s s β-1 ds β (s -s β ) γ β+1 2 exp    i β-1 b=0 (p -q + Φ b )[θ(s b+1 ) -θ(s b )]    β-1 b=0 (s b+1 -s b ) γ b+1 2 e -(p-q+Φ b ) 2 K θ (s b+1 -s b ) 2 ⌊ γ b+1 2 ⌋ a=0 γ b+1 ! i(p -q + Φ b ) K θ (s b+1 -s b ) γ b+1 -2a a!(γ b+1 -2a)!2 a     
where the sum over γ, i.e. γ , is the sum over all the vectors γ =

[γ 1 , γ 2 , • • • , γ β , γ β+1 ],
where γ b (b = 1, • • • , β) are positive integers and γ β+1 is a positive integer with even parity. Additionally, they satisfy the constraint: β+1 b=1 γ b = r. The expressions provided in the previous section allow us to compute any order order statistics for the considered Brownian motion. These expressions contain multiple integrals on the active trajectory. In some cases, these integrals can be easily solved. This is certainly the case when the expression of θ(s) is a linear function of the curve length s. Let us refer to this special case.

The case of circular active trajectories

First of all, we remark that, when the orientation linearly depends on the curve length, the particle accomplishes a circular trajectory. From the third equation in (3), by only considering the deterministic part, we obtain that θ(s) is linear in s when µ(s) = µ 0 , namely, when the ratio of the angular and linear speed is constant. We have:

θ(s) = θ 0 + µ 0 s (19)
and, the corresponding active trajectory, is a circumference with radius 1 µ 0 (note that when the angular speed vanishes, the radius becomes infinite and the active trajectory becomes a straight line). The computation of the integrals in ( 17) and ( 18) is immediate.

We conclude by providing the expression of the following mean values: x(s) , y(s) and D(s) 2 . Let us introduce the following complex quantity:

z ≡ - K θ 2 + iµ 0 (20)
From ( 8) and ( 9) we obtain x(s) + i y(s) = e iθ 0 s 0 ds ′ e zs ′ = e iθ 0 e zs -1 z . Hence:

x(s) = cos θ 0 ℜ {f (z, s)} -sin θ 0 ℑ {f (z, s)} (21)
y(s) = cos θ 0 ℑ {f (z, s)} + sin θ 0 ℜ {f (z, s)} (22)
where f (z, s) ≡ e zs -1 z and the symbols ℜ{•} and ℑ{•} are the real and the imaginary part of a given complex quantity. Finally, from [START_REF] Oksendal | Stochastic Differential Equations: An Introduction with Applications[END_REF] we obtain:

D(s) 2 = K r s + 2ℜ f (z, s) -s z (23)
Note that this last mean value is independent of the initial orientation, as expected.

Obviously, this property holds for any active trajectory and can directly be remarked in [START_REF] Oksendal | Stochastic Differential Equations: An Introduction with Applications[END_REF].

Conclusion

In this paper we derived the statistics, up to any order, for an overdamped 2D-Brownian motion of an active and non-holonomic particle. The chosen kinematic constraint is modelled by the unicycle differential equation and it is a very general constraint for a 2D motion. According to this model, the particle can freely rotates but only one direction for the shift is allowed. The considered active motion is very general. The expressions here provided for the statistics hold for any active trajectory that satisfies the mentioned kinematic constraint. The expressions contain multiple integrals over the active trajectory. In the case when the ratio between the angular and linear speed is constant, these integrals can be easily computed. In particular, some expressions for the statistics up to the second order are provided. We want to remark that this paper provides the analytical expressions for the statistics of a not trivial Brownian motion up to any order. To the best of our knowledge, this has never been done in the past. Deriving any-order moment corresponds to analytically derive the probability density distribution [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF]. In this sense, the paper provides the analytical solution to a not trivial Fokker-Plank equation, for which the detailed balance does not hold.

Appendix A. Computation of u(s) p v(s) q

We have:

u(s) p v(s) q = lim N →∞ j 1 •••jpk 1 •••kq e i(θ j 1 +•••+θ jp -θ k 1 -•••-θ kq ) (A.1) (δs + ǫ j 1 ) • • • (δs + ǫ jp )(δs + ǫ k 1 ) • • • (δs + ǫ kq )
where each index goes from 1 to N . Let us focus our attention on the quantity (δs + ǫ j 1 )

• • • (δs + ǫ jp )(δs + ǫ k 1 ) • • • (δs + ǫ kq ) .
We can write this quantity as the sum of p + q + 1 terms, i.e., δs p+q C 0 + δs p+q-1 C 1 + δs p+q-2 C 2 + • • • + δsC p+q-1 + C p+q . We trivially have C 0 = 1. Concerning the remaining coefficients, we remark, first of all, that only the ones with even index are different from 0. In other words,

C 2n+1 = 0, n = 0, 1, • • • , ⌊ p+q 2 ⌋
(where ⌊x⌋ is the largest integer not greater than x). Let us compute C 2n . This coefficient is the sum of p+q 2n elements, each of them being the average of a product of 2n terms "ǫ". On the other hand, since the exponential in (A.1) is symmetric with respect to the change j a ↔ j a ′ , a, a ′ = 1, • • • , p and with respect to the change

k b ↔ k b ′ , b, b ′ = 1, • • • , q but it
is not symmetric with respect to the change j a ↔ k b , we need to separate the p+q 2n elements in several groups. Specifically, let us denote by l the number of ǫ with index of type k. The number of elements belonging to this group is p 2n-l q l . Note that we set a b = 0 when b > a. Note also that 2n l=0 p 2n-l q l = p+q 2n . We now remark that, because of the statistical properties of ǫ, the average of the product of 2n terms ǫ is different from zero only when their indexes are equal two-by-two. Hence, we have to consider all the combinations of products of terms ǫ, whose indexes are equal two-by-two and that differ for at least one pair of indexes. On the other hand, for a given element in the group characterized by a given n and l, the effect in (A.1) depends on the number of pairs that are hetero (i.e., with an index of type j and one of type k). Let us denote by m the number of pairs that are hetero. Note that m has the same parity of l. Indeed, by definition we have l indexes of type k. Additionally, we are using m indexes of type j and m indexes of type k to make m hetero pairs. Hence, lm indexes of type k remain and, with them, we have to make l-m 2 homo pairs of type k. Similarly, we have to make 2n-l-m 2 homo pairs of type j. Hence, l-m 2 must be integer. We compute the number of combinations of the products of 2n terms ǫ, with l indexes of type k, where the indexes are equal two-by-two, with m pairs that are hetero and that differ for at least one pair. We obtain: 2n-l m l m m!(2n-l-m-1)!!(l-m-1)!!. Following this, we can write (A.1) as follows:

u(s) p v(s) q = lim N →∞ {j}p{k}q ⌊ p+q 2 ⌋ n=0 δs p+q-2n 2n l=0 p 2n -l q l (A.2) • min(l,2n-l) m=0 odd/even 2n -l m l m m!(2n -l -m -1)!!(l -m -1)!! ǫ 2 j 1 δ j 1 j 2 ǫ 2 j 3 δ j 3 j 4 • • • ǫ 2 j 2n-l-m-1 δ j 2n-l-m-1 j 2n-l-m ǫ 2 k 1 δ k 1 k 2 ǫ 2 k 3 δ k 3 k 4 • • • ǫ 2 k l-m-1 δ k l-m-1 k l-m ǫ 2 k l-m+1 δ k l-m+1 j 2n-l-m+1 • • • ǫ 2 k l δ k l j 2n-l e i(θ j 1 +•••+θ jp -θ k 1 -•••-θ kq )
where, for the brevity sake, we denoted by {j}p{k}q the sum

N j 1 •••jpk 1 •••kq=1
. Each average ǫ 2 provides K r δs. Hence, all together, they provide K n r δs n . Additionally, the number of Kronecker deltas is n. Hence, in the limit of N → ∞ for each value of n we get a multiple integral of dimension p + qn. Note that, when the indexes are equal h-by-h (h > 2), the result in the limit N → ∞ vanishes since, the power of δs, is larger than the number of sums. By a direct computation in (A.2) we obtain:

u(s) p v(s) q = ⌊ p+q 2 ⌋ n=0 K n r 2n l=0 p 2n -l q l min(l,2n-l) m=0 odd/even 2n -l m l m (A.3) m!(2n -l -m -1)!!(l -m -1)!! s m lim N →∞ {j s }σ{j d }ρ{k s }χ{k d }η δs β exp i[θ j s 1 + • • • + θ j s σ + 2(θ j d 1 + • • • + θ j d ρ ) -(θ k s 1 + • • • + θ k s χ ) -2(θ k d 1 + • • • + θ k d η )] where ρ ≡ 2n-l-m 2 , σ ≡ p-(2n-l), η ≡ l-m
2 , χ ≡ q-l and β ≡ ρ+σ+η+χ = p+q-n-m. We must compute the average of the exponential in (A.3) and then we compute the limit N → ∞. We remark that the various θ in the exponential contain random quantities (i.e., the δθ at different time steps). In order to proceed we have to separate all the random quantities that are independent. We start this separation by redefining the indexes in the sum {j s }σ{j d }ρ{k s }χ{k d }η . Specifically, we consider the new indexes j s j d k s k d that differ from j s j d k s k d since they are ordered (in increasing order). For instance, j

s 1 < j s 2 < • • • < j s σ .
Hence, the last sum in (A.3) can be replaced with {j s }σ{j d }ρ{k s }χ{k d }η → ρ!σ!η!χ! {j s }σ{j d }ρ{k s }χ{k d }η . The four types of indexes are not ordered among them. Hence, the sum includes all the possible combinations that maintain the order only restricted to a single index type. For instance, a possible combination is: that Φ β = -α. According to this, the last sum in (A.3) can be written as follows:

k d 1 < k d 2 < j s 1 < k s 1 < j s 2 < k d 3 < j d 1 < j d 2 < • • •. We introduce the β ordered indexes w 1 < w 2 < • • • < w β .
e iαθ 0 ρ!σ!η!χ! c {w} β δs β exp i β-1 b=0 (α + Φ b )(∆ w b+1 w b + ∆ w b+1
w b ) . Hence, we have:

u(s) p v(s) q = ⌊ p+q 2 ⌋ n=0 K n r 2n l=0 p 2n -l q l min(l,2n-l) m=0 odd/even 2n -l m l m (A.4) m!(2n -l -m -1)!!(l -m -1)!! s m e iαθ 0 ρ!σ!η!χ! c lim N →∞ {w} β δs β exp    i β-1 b=0 (α + Φ b )(∆ w b+1 w b + ∆ w b+1 w b )   
Now we can compute the average since we were able to separate all the independent quantities. We have:

exp    i β-1 b=0 (α + Φ b )(∆ w b+1 w b + ∆ w b+1 w b )    = (A.5) = exp    i β-1 b=0 (α + Φ b )∆ w b+1 w b    β-1 b=0 exp i (α + Φ b )∆ w b+1 w b = = exp    i β-1 b=0 (α + Φ b )∆ w b+1 w b    β-1 b=0 exp - (α + Φ b ) 2 (w b+1 -w b )K θ δs 2 = = exp    β-1 b=0 i(α + Φ b )∆ w b+1 w b - (α + Φ b ) 2 (w b+1 -w b )K θ δs 2   
We used the equality (η = N (0, 1)):

e Aη = e A 2 2
(A.6)

with A = i(α + Φ b )σ b+1 and σ 2 b+1 = (w b+1w b )K θ δs. By substituting equation (A.5) in (A.4) and by taking the limit (N → ∞), we finally obtain:

u(s) p v(s) q = ⌊ p+q 2 ⌋ n=0 K n r 2n l=0 p 2n -l q l min(l,2n-l) m=0 odd/even 2n -l m l m (A.7) m!(2n -l -m -1)!!(l -m -1)!! s m e iαθ 0 ρ!σ!η!χ! c s 0 ds 1 s s 1 ds 2 • • • s s β-1 ds β exp    β-1 b=0 i(α + Φ b )[θ(s b+1 ) -θ(s b )] - (α + Φ b ) 2 (s b+1 -s b )K θ 2   
which coincides with [START_REF] Peter | Brownian motors: noisy transport far from equilibrium[END_REF].

Appendix B. Computation of u(s) p v(s) q θr This computation follows the same initial steps carried out in appendix Appendix A. We obtain an expression equal to the one given in (A.4) but the mean value at the end must be replaced with exp i β-1 b=0 (α + Φ b )(∆

w b+1 w b + ∆ w b+1 w b ) • θr .
The deterministic part can be factorized out the mean value. We need to calculate: which coincides with [START_REF] Risken | The Fokker-Planck Equation Methods of Solution and Applications[END_REF].

exp    i β-1 b=0 (α + Φ b )∆ w b+1 w b    • θr = θr β-1 b=0 exp i(α + Φ b )∆

Figure 1 .

 1 Figure 1. ...

Figure 2 .

 2 Figure 2. ...

+ ∆ w 1 0 2 w 1 ) 2 w 1 ) 2 w 1 ) 2 w 1 ) 2 w 1 )

 12121212121 Additionally, let us denote with ∆ b a ≡ b c=a+1 δθ c and ∆ b a ≡ b c=a+1 δθ c = θ(bδs)θ(aδs), where θ(s) is the active trajectory. Finally, let us define α ≡ 2ρ + σχ -2η. The sum in the exponential in (A.3) contains α(∆ w 1 0). Then, depending on the combination of the indexes js j d k s k d , we have a different result for the (∆ w 2 w 1 + ∆ w . Specifically, if w 1 = j d 1 the sum in the exponential contains (α -2)(∆ w 2 w 1 + ∆ w . If w 1 = j s 1 it contains (α -1)(∆ w 2 w 1 + ∆ w . If w 1 = k s 1 it contains (α + 1)(∆ w 2 w 1 + ∆ w . If w 1 = k d 1 it contains (α + 2)(∆ w 2 w 1 + ∆ w . We introduce the vector c ≡ [c 1 , • • • , c β ] T whose entries are -2, -1,1 and 2. Specifically, it contains η entries equal to 2, χ equal to 1, σ equal to -1 and ρ equal to -2. It is easy to realize that we have β χ β-χ σ ρ+η η vectors c. Finally, we define Φ b ≡ b a=1 c a . Note

γγ(s b+1 -s b ) γ b+1 2 eγ(s b+1 -s b ) γ b+1 2 eγ

 22 where w β+1 ≡ N and the sum over γ is the sum over all the vectors γ = [γ 1 , • • • , γ β+1 ], where γ b (b = 1, • • • , β + 1) are positive integers satisfying the constraint: β+1 b=1 γ b = r. Hence, we obtain θr We use the following two standard results for a normal distribution: β+1 -1)!! γ β+1 even (B.2) andexp i(α + Φ b )∆ b+1 ! (i(α + Φ b )σ b+1 ) γ b+1 -2a a!(γ b+1 -2a)!2 a with σ 2 b+1 ≡ K θ δs(w b+1w b ). Equation (B.3) is obtained starting from (A.6) and by differentiating γ b+1 times with respect to A. Hence, the expression in (B.1), becomes:θr β-1 b=0 exp i(α + Φ b )∆ w b+1 w b = γ r!σ γ β+1 β+1 (γ β+1 -1)!! β b=0 γ b+1 ! b+1 ! (i(α + Φ b )σ b+1 ) γ b+1 -2aa!(γ b+1 -2a)!2 a sum over γ only includes the vectors γ whose last entry is even. By taking the limit N → ∞, this expression becomes:-(α+Φ b ) 2 K θ (s b+1 -s b ) b+1 ! i(α + Φ b ) K θ (s b+1s b ) γ b+1 -2aa!(γ b+1 -2a)!2 a obtain:u(s) p v(s) q θ(s) r = 2nlm -1)!!(lm -1)!! s m e iαθ 0 ρ!σ!η!χ! Φ b )[θ(s b+1 )θ(s b )] -(α+Φ b ) 2 K θ (s b+1 -s b ) b+1 ! i(α + Φ b ) K θ (s b+1s b ) γ b+1 -2a a!(γ b+1 -2a)!2 a

‡ Note that, deriving any-order moment, corresponds to deriving the probability density distribution[START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF].