# Far-From-Equilibrium Time Evolution between Two Gamma Distributions

Eun-Jin Kim, Lucille-Marie Tenkès, Rainer Hollerbach, Ovidiu Radulescu

► To cite this version:

Eun-Jin Kim, Lucille-Marie Tenkès, Rainer Hollerbach, Ovidiu Radulescu. Far-From-Equilibrium Time Evolution between Two Gamma Distributions. Entropy, 2017, 19 (10), pp.511. 10.3390/e19100511. [https://hal.science/hal-03644310](https://hal.science/hal-03644310)

HAL Id: hal-03644310

https://hal.science/hal-03644310

Submitted on 27 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Far-From-Equilibrium Time Evolution between Two Gamma Distributions

Eun-jin Kim<fnref n="1" />,*, Lucille-Marie Tenkès<fnref n="1" />,<fnref n="2" />, Rainer Hollerbach<fnref n="3" /> and Ovidiu Radulescu<fnref n="4" />

Received: 18 August 2017; Accepted: 21 September 2017; Published: 22 September 2017

**Abstract:** Many systems in nature and laboratories are far from equilibrium and exhibit significant fluctuations, invalidating the key assumptions of small fluctuations and short memory time in or near equilibrium. A full knowledge of Probability Distribution Functions (PDFs), especially time-dependent PDFs, becomes essential in understanding far-from-equilibrium processes. We consider a stochastic logistic model with multiplicative noise, which has gamma distributions as stationary PDFs. We numerically solve the transient relaxation problem and show that as the strength of the stochastic noise increases, the time-dependent PDFs increasingly deviate from gamma distributions. For sufficiently strong noise, a transition occurs whereby the PDF never reaches a stationary state, but instead, forms a peak that becomes ever more narrowly concentrated at the origin. The addition of an arbitrarily small amount of additive noise regularizes these solutions and re-establishes the existence of stationary solutions. In addition to diagnostic quantities such as mean value, standard deviation, skewness and kurtosis, the transitions between different solutions are analysed in terms of entropy and information length, the total number of statistically-distinguishable states that a system passes through in time.

**Keywords:** non-equilibrium statistical mechanics; gamma distribution; stochastic processes; Fokker-Planck equation; fluctuations and noise

## 1. Introduction

In classical statistical mechanics, the Gaussian (or normal) distribution and mean-field type theories based on such distributions have been widely used to describe equilibrium or near equilibrium phenomena. The ubiquity of the Gaussian distribution stems from the central limit theorem that random variables governed by different distributions tend to follow the Gaussian distribution in the limit of a large sample size [1–3]. In such a limit, fluctuations are small and have a short correlation time, and mean values and variance completely describe all different moments, greatly facilitating analysis.

Many systems in nature and laboratories are however far from equilibrium, exhibiting significant fluctuations. Examples are found not only in turbulence in astrophysical and laboratory plasmas, but also in forest fires, the stock market and biological ecosystems [4–23]. Specifically, anomalous (much larger than average values) transport associated with large fluctuations in fusion plasmas can degrade the confinement, potentially even terminating fusion operation [6]. Tornadoes are rare, large amplitude events, but can cause very substantial damage when they do occur. In biology, the pioneering work of Delbrück on bacteriophages showed that viruses replicate in strongly fluctuating bursts [24]. The fluctuations of the burst amplitudes were explained mathematically by stochastic

<footnote n="1">¹ School of Mathematics and Statistics, University of Sheffield, Sheffield S3 7RH, UK; lucillemarie.tenkes@gmail.com</footnote>

<footnote n="2">² ENSTA ParisTech Université Paris-Saclay, 828 boulevard des Maréchaux, 91120 Palaiseau, France</footnote>

<footnote n="3">³ Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK; R.Hollerbach@leeds.ac.uk</footnote>

<footnote n="4">⁴ DIMNP-UMR 5235 CNRS, Université de Montpellier, Place Eugène Bataillon, 34095 Montpellier, France; ovidiu.radulescu@umontpellier.fr</footnote>

<footnote>* Correspondence: e.kim@shef.ac.uk; Tel.: +44-114-222-3876</footnote>
autocatalytic reaction models first introduced in [25]. Delbrück's autocatalytic models predict discrete negative-binomial distributions, which can be well approximated by gamma distributions when the average number of particles is large. Furthermore, gene expression and protein productions, which used to be thought of as smooth processes, have also been observed to occur in bursts, leading to negative binomial and gamma distributed protein copy numbers (e.g., [19–23]). Such rare events of large amplitude (called intermittency) can dominate the entire transport even if they occur infrequently [8,26]. They thus invalidate the assumption of small fluctuations with short correlation time, making mean value and variances meaningless. Therefore, to understand the dynamics of a system far from equilibrium, it is crucial to have a full knowledge of Probability Distribution Functions (PDFs), including time-dependent PDFs [27].

The consequences of strong fluctuations in far-from-equilibrium systems are multiple. In physics, far-from-equilibrium fluctuations produce dissipative patterns, shift or wipe out phase transitions, etc. In economics, finance and actuarial science, strong fluctuations are important issues of risk evaluation. In biology, strong fluctuations generate phenotypic heterogeneity that helps multicellular organisms or microbial populations to adapt to changes of the environment by so-called "bet-hedging" strategies. In such a strategy, only a part of the cell population adapts upon emergence of new environmental conditions. The remaining part retains the memory of the old conditions and is thus already adapted if environmental conditions revert to previous ones [28]. Exceptional behaviour can also rescue cell subpopulations from drug-induced lethal conditions, thus generating drug resistance [29]. In particular, because of the skewness and exponential tail of the gamma distribution, gamma-distributed populations contain a significant proportion of individuals with an exceptionally high phenotypic variation. Bet-hedging being a dynamic phenomenon, it is important, for biological studies, to be able to predict not only steady-state, but also time-dependent distributions.

Obtaining good quality PDFs is often very challenging, as it requires a sufficiently large number of simulations or observations. Therefore, a PDF is usually constructed by averaging data from a long time series and is thus stationary (independent of time). Unfortunately, such stationary PDFs miss crucial information about the dynamics/evolution of non-equilibrium processes (e.g., tumour evolution). Theoretical prediction of time-dependent PDFs has proven to be no less challenging due to the limitation in our understanding of nonlinear stochastic dynamical systems, as well as the complexity in the required analysis.

Spectral analysis, for example, using theoretical tools similar to those used in quantum mechanics (e.g., raising and lower operators) is useful (e.g., [1]), but the summation of all eigenfunctions is necessary for time-dependent PDFs far from equilibrium. Various different methodologies have also been developed to obtain approximate PDFs, such as the variational principle, the rate equation method or the moment method [30–35]. In particular, the rate equation method [31,32] assumes that the form of a time-dependent PDF during the relaxation is similar to that of the stationary PDF and thus approximates a time-dependent PDF during transient relaxation by a PDF having the same functional form as a stationary PDF, but with time-varying parameters.

In this work, we show that this assumption is not always appropriate. We consider a stochastic logistic model with multiplicative noise. We show that for fixed parameter values, the stationary PDFs are always gamma distributions (e.g., [36,37]), one of the most popular distributions used in fitting experimental data. However, we find numerically that the time-dependent PDFs in transitioning from one set of parameter values to another are significantly different from gamma distributions, especially for strong stochastic noise. For sufficiently strong multiplicative noise, it is necessary to introduce additive noise, as well, to obtain stationary distributions at all. We note that in inferential statistics, gamma distributions facilitate Bayesian model learning from data, as a gamma distribution is a conjugate prior to many likelihood functions. It is therefore interesting to test whether models with stationary gamma distributions also have time-dependent gamma distributions.
## 2. Stochastic Logistic Model

We consider the logistic growth with a multiplicative noise given by the following Langevin equation:

$$ \frac{dx}{dt} = (\gamma + \zeta)x - \epsilon x^2, \qquad (1) $$

where $x$ is a random variable and $\zeta$ is a stochastic forcing, which for simplicity, can be taken as a short-correlated random forcing as follows:

$$ \langle \zeta(t) \zeta(t') \rangle = 2D\delta(t - t'). \qquad (2) $$

In Equation (2), the angular brackets represent the average over $\zeta$, $\langle \zeta \rangle = 0$, and $D$ is the strength of the forcing. $\gamma$ is the control parameter in the positive feedback, representing the growth rate of $x$, while $\epsilon$ represents the efficiency in self-regulation by a negative feedback. $\gamma x - \epsilon x^2$ can be considered as the gradient of the potential $V$ as $\gamma x - \epsilon x^2 = -\frac{\partial V}{\partial x}$, where $V = -\frac{\gamma}{2}x^2 + \frac{\epsilon}{3}x^3$. Thus, $V$ has its minimum value at $x = \frac{\gamma}{\epsilon}$. When $\zeta = 0$, $x = \frac{\gamma}{\epsilon}$ (the carrying capacity) is a stable equilibrium point since $\partial_{xx}V|_{x=\gamma/\epsilon} = \gamma > 0$; $x = 0$ is an unstable equilibrium point since $\partial_{xx}V|_{x=0} = -\gamma < 0$.

The multiplicative noise in Equation (1) shows that the linear growth rate contains the stochastic part $\zeta$. This model is entirely phenomenological, and $x$ can be interpreted as the size of a critical physical phenomenon (vortex, tornado, etc.), stock market, number of biological cells, viruses and proteins. It is reminiscent of Delbrück's autocatalytic processes [25], but is different from these in many ways, the most important being the lack of discreteness and the possibility of reaching a steady-state due to the finite capacity of logistic growth. We will show in the following that in spite of these differences, our model is capable of producing large fluctuations.

By using the Stratonovich calculus [2,3,38], we can obtain the following Fokker-Planck equation for the PDF $p(x, t)$ (see Appendix A for details):

$$ \frac{\partial}{\partial t} p(x, t) = -\frac{\partial}{\partial x} \left[ (\gamma x - \epsilon x^2) p(x, t) \right] + D \frac{\partial}{\partial x} \left[ x \frac{\partial}{\partial x} \left[ x p(x, t) \right] \right] \qquad (3) $$

corresponding to the Langevin Equation (1). By setting $\partial_t p = 0$, we can analytically solve for the stationary PDFs as:

$$ p(x) = \frac{b^a}{\Gamma(a)} x^{a-1} e^{-bx}, \qquad (4) $$

which is the well-known gamma distribution. The two parameters $a$ and $b$ are given by $a = \gamma/D$ and $b = \epsilon/D$. The mean value and variance of the gamma distribution are found to be:

$$ \langle x \rangle = \frac{a}{b} = \frac{\gamma}{\epsilon}, \qquad \text{Var}(x) = \sigma^2 = \langle (x - \langle x \rangle)^2 \rangle = \frac{a}{b^2} = \frac{\gamma D}{\epsilon^2}, \qquad (5) $$

where $\sigma = \sqrt{\text{Var}(x)}$ is the standard deviation. We recognise $\langle x \rangle$ as the carrying capacity for a deterministic system with $\zeta = 0$. It is useful to note that $\langle x \rangle$ is given by the linear growth rate scaled by $\epsilon$, while $\text{Var}(x)$ is given by the product of the linear growth rate and the diffusion coefficient, each scaled by $\epsilon$. That is, the effect of stochasticity should be measured relative to the linear growth rate.

Therefore, the case of small fluctuations is modelled by values of $D$ small compared with $\gamma$ and $\epsilon$. In such a limit, $a$ and $b$ are large, making $\sqrt{\text{Var}(x)} \ll \langle x \rangle$ in Equation (5). That is, the width of the PDF is much smaller than its mean value. In this limit, Equation (4) reduces to a Gaussian distribution. To show this, we express Equation (4) in the following form:

$$ p \equiv \frac{b^a}{\Gamma(a)} e^{-f(x)}, \qquad (6) $$
where $f(x) = bx - (a - 1) \ln x$. For large $b$, we expand $f(x)$ around the stationary point $x = x_p$ where $\partial_x f(x) = 0 = b - (a - 1)/x$ up to the second order in $x - x_p$ to find:

$$x_p = \frac{a-1}{b} \sim \frac{a}{b'}, \qquad f(x = x_p) \sim a \left(1 - \ln \frac{a}{b}\right), \qquad (7)$$

$$f(x) \sim f(x_p) + \frac{1}{2}(x - x_p)^2 \partial_{xx} f(x) \Big|_{x=x_p} = a \left(1 - \ln \frac{a}{b}\right) + \frac{b^2}{2a} \left(x - \frac{a}{b}\right)^2. \qquad (8)$$

Here, $a \gg 1$ was used. Using Equation (8) in Equation (6) then gives us:

$$p \propto \exp \left[ -\frac{b^2}{2a} \left( x - \frac{a}{b} \right)^2 \right] \propto \exp \left[ -\beta (x - \langle x \rangle)^2 \right], \qquad (9)$$

which is a Gaussian PDF with mean value $\langle x \rangle$. Here, $\beta = 1/\text{Var}(x)$ is the inverse temperature, and the variance $\text{Var}(x)$ is given by Equation (5). Therefore, for a sufficiently small $D$, the gamma distribution is approximated as a Gaussian PDF, which is consistent with the central limit theorem as small $D$ corresponds to small fluctuations and large system size. See also [39] for a different derivation.

As $D$ increases, the Gaussian approximation becomes increasingly less valid. Indeed, even the gamma distribution becomes invalid asymptotically, when $t \to \infty$, if $D > \gamma$; according to Equation (4), having $a < 1$ yields $\lim_{x \to 0} p = \lim_{x \to 0} \frac{\partial p}{\partial x} = \infty$. However, from the full time-dependent Fokker-Planck Equation (3), one finds that if the initial condition satisfies $p = 0$ at $x = 0$, then $p(x = 0)$ will remain zero for all later times. As we will see, the resolution to this seeming paradox is that no stationary distribution is ever reached for $D > \gamma$, but instead, the peak approaches ever closer to $x = 0$, without ever reaching it.

If we are interested in obtaining stationary solutions even when $D > \gamma$, one way to achieve that is to return to the original Langevin Equation (1), but now include a further additive stochastic noise $\eta$:

$$\frac{dx}{dt} = (\gamma + \zeta)x - \epsilon x^2 + \eta, \qquad (10)$$

where $\zeta$ and $\eta$ are uncorrelated and $\eta$ satisfies $\langle \eta(t)\eta(t') \rangle = 2Q\delta(t - t')$. The new version of the Fokker-Planck Equation (3) then becomes:

$$\frac{\partial}{\partial t} p = -\frac{\partial}{\partial x} \left[ (\gamma x - \epsilon x^2) p \right] + D \frac{\partial}{\partial x} \left[ x \frac{\partial}{\partial x} [x p] \right] + Q \frac{\partial^2}{\partial x^2} p, \qquad (11)$$

which has stationary solutions given by:

$$\ln p(x) = \int \frac{(\gamma - D)x - \epsilon x^2}{Dx^2 + Q} dx. \qquad (12)$$

This integral can be evaluated analytically, but the final form is not particularly illuminating. The only point to note is that for non-zero $Q$, the denominator is never zero, even for $x \to 0$, which avoids any possible singularities at the origin. For $\gamma > D$ and $Q \ll D$, the solutions are also essentially indistinguishable from the previous gamma distribution (4). The only significant effect of including $\eta$ therefore is to avoid the previous difficulties at the origin when $D > \gamma$.

As we have seen, both Fokker-Planck Equations (3) and (11) can be solved exactly for their stationary solutions. This is unfortunately not the case regarding time-dependent solutions, where no closed-form analytic solutions exist (see Appendix B for the extent to which analytic progress can be made). We therefore developed finite-difference codes, second-order accurate in both space and time. Most aspects of the numerics are standard and similar to previous work [40–42]. The only point that requires discussion are the boundary conditions. As noted above, for (3), the equation itself states that $p = 0$ at $x = 0$ is the appropriate boundary condition, provided only that the initial condition also
satisfies this. In contrast, for (11), the appropriate boundary condition is $\frac{\partial}{\partial x}p = 0$ at $x = 0$. To derive this boundary condition for (11), we simply integrate (11) over the range $x = [0, \infty]$ and require that the total probability should always remain one, so that $\frac{d}{dt} \int p dx = 0$. Regarding the outer boundary, choosing some moderately large outer value for $x$ and then imposing $p = 0$ there was sufficient. Resolutions up to $10^6$ grid points were used, and results were carefully checked to ensure they were independent of the grid size, time step and precise choice of outer boundary.

### 3. Diagnostics

Once the time-dependent solutions are computed, we can analyse them using a number of diagnostics. First, we can evaluate the mean value $\langle x \rangle$ and standard deviation $\sigma$ from (5). Next, to explore the extent to which the time-dependent PDFs differ from gamma distributions, we can simply compare them with 'equivalent' gamma distributions and compute the difference. That is, given $\langle x \rangle$ and $\sigma$, the gamma distribution $p_{\text{equiv}}$ having the same mean and variance would have as its two parameters $a = \langle x \rangle^2 / \sigma^2$ and $b = \langle x \rangle / \sigma^2$. With these values, we define:

$$ \text{Difference} = \int |p - p_{\text{equiv}}| dx \qquad (13) $$

to measure how different the actual time-dependent PDF is from its equivalent gamma distribution.

Two other familiar quantities often useful in analysing PDFs are the skewness and kurtosis, defined by:

$$ \text{Skewness} = \frac{\langle (x - \langle x \rangle)^3 \rangle}{\sigma^3}, \qquad \text{Kurtosis} = \frac{\langle (x - \langle x \rangle)^4 \rangle}{\sigma^4} - 3. \qquad (14) $$

Skewness measures the extent to which a PDF is asymmetric about its peak, whereas kurtosis measures how concentrated a PDF is in the peak versus the tails, relative to a Gaussian having the same variance (the -3 is included in the definition of the kurtosis to ensure that a Gaussian would yield zero). For gamma distributions, one finds analytically that the skewness is $2\sqrt{D/\gamma}$, and the kurtosis is $6D/\gamma$. Comparing the skewness and kurtosis of the time-dependent PDFs with these formulas is therefore another useful way of quantifying how similar or different they are from gamma distributions.

Another quantity that can be useful is the so-called differential entropy as a measure of order versus disorder (as entropy always is):

$$ S = - \int p \ln p dx. \qquad (15) $$

In particular, we expect $S$ to be small for localised PDFs and large for spread out ones (e.g., [40–43]). For unimodal PDFs as the ones studied here, entropy and standard deviation are typically comparably good measures of localization, but for bimodal peaks, entropy can be significantly better [42]. For the gamma distribution in Equation (4), the differential entropy can be shown to be given by:

$$ S = a - \ln b + \ln(\Gamma(a)) + (1 - a)\psi(a), \qquad (16) $$

where $\psi(a) = \frac{d \ln(\Gamma(a))}{dx}|_{x=a}$ is the double gamma function.

Our final diagnostic quantity is what is known as information length. Unlike all the previous diagnostics, which are simply evaluated at any instant in time, but otherwise do not involve $t$, information length is a Lagrangian quantity, explicitly concerned with the full time-history of the evolution of a given PDF. It is thus ideally suited to understanding time-dependent PDFs. Very briefly, we begin by defining:

$$ \mathcal{E} \equiv \frac{1}{[\tau(t)]^2} = \int \frac{1}{p(x,t)} \left[ \frac{\partial p(x,t)}{\partial t} \right]^2 dx. \qquad (17) $$

Note how $\tau$ has units of time and quantifies the correlation time over which the PDF changes, thereby serving as a time unit in statistical space. Alternatively, $1/\tau$ quantifies the (average) rate of change of
information in time. $\mathcal{E}$ is due to the change in either width (variance) of the PDF or the mean value, which are determined by $\gamma$, $D$ and $\epsilon$ for the gamma distribution (e.g., see Equation (4)). In the standard Brownian motion, the mean value is zero so that $\mathcal{E}$ is due to the change in the variance of a PDF.

The total change in information between initial and final times, zero and $t$ respectively, is then defined by measuring the total elapsed time in units of $\tau$ as:

$$ \mathcal{L}(t) = \int_{0}^{t} \frac{dt_1}{\tau(t_1)} = \int_{0}^{t} \sqrt{\int dx \frac{1}{p(x, t_1)} \left[ \frac{\partial p(x, t_1)}{\partial t_1} \right]^2} dt_1. \quad (18) $$

This information length $\mathcal{L}$ measures the total number of statistically distinguishable states that a system evolves through, thereby establishing a distance between the initial and final PDFs in the statistical space. Note that $\mathcal{L}$ by construction is a continuous variable and thus measures the total “number” of statistically-different states as a continuous number. See also [40–48] for further applications and theoretical background of $\mathcal{E}$ and $\mathcal{L}$.

# 4. Results

## 4.1. $\gamma > D$

We start with the case $\gamma > D$, where Equation (3) has stationary solutions, given by (4). Keeping $\epsilon$ and $D$ fixed, we then switch $\gamma$ back and forth between two values, in the following sense: Take the gamma distribution (4) corresponding to one value; call it $\gamma_1$; and use that as the initial condition to solve (3) with the other value; call it $\gamma_2$. We then interchange $\gamma_1$ and $\gamma_2$ to complete the pair of “inward” and “outward” processes. Such a pair can be thought of as an order/disorder phase transition [40,41], caused for example by cyclically adjusting temperature in an experiment. This protocol is also inspired from adaptation of a biological system. During adaptation, a model parameter can be abruptly changed in response to the change of environmental conditions, for instance a particle replication parameter $\gamma$, but the resulting changes can be extremely heterogeneous in the population.

Figure 1 shows the result of switching $\gamma$ between $\gamma_1 = 0.5$ and $\gamma_2 = 0.05$, for fixed $\epsilon = 1$ and $D = 0.02$ (one of the three parameters $\epsilon$, $\gamma$ and $D$ can of course always be kept fixed by rescaling the entire equation, so throughout this entire section, we keep $\epsilon = 1$ fixed and focus on how the various quantities depend on $\gamma$ and $D$). We immediately see that the inward and outward processes behave differently. When $\gamma$ is decreased and the peak therefore moves inward, the PDF is relatively narrow, and the peak amplitude is monotonically increasing. When $\gamma$ is switched from 0.05 back to 0.5, the PDF is much broader, and the peak amplitude in the intermediate stages is less than either the initial or final gamma distributions.

Figure 2 shows how $\langle x \rangle$, $\mathcal{E}$ and $\mathcal{L}$ vary as functions of time, for the three values $D = 0.01, 0.02, 0.04$. For $\langle x \rangle$, the movement from 0.5 to 0.05 is somewhat slower than the reverse process; but, both processes occur on a similar time-scale, and both are essentially independent of $D$. This is in contrast with other Fokker-Planck systems where the magnitude of the diffusion coefficient can have a very strong influence on the equilibration time-scales [40,41].
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) shows the result of switching $\gamma = 0.5 \rightarrow 0.05$; (b) $\gamma = 0.05 \rightarrow 0.5$, both at fixed $\epsilon = 1$ and $D = 0.02$. The initial (red) and final (blue) gamma distributions are shown as heavy lines. The four intermediate lines are when the time-dependent solutions have $\langle x \rangle = 0.1, 0.2, 0.3, 0.4$. The arrows are a reminder of the direction of motion, inward on the left and outward on the right.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) shows $\langle x \rangle$ as a function of time; (b) shows $\mathcal{E} \cdot D$ (to indicate the $\mathcal{E} \sim D^{-1}$ scaling); (c) shows $\mathcal{L} \cdot D^{1/2}$ (to indicate the $\mathcal{L} \sim D^{-1/2}$ scaling). Solid lines denote $\gamma = 0.5 \rightarrow 0.05$, dashed lines indicate the reverse. Each solid or dashed “line” is in fact three, occasionally barely distinguishable, lines with $D = 0.01, 0.02, 0.04$. The dots on the $\langle x \rangle$ curves correspond to the PDFs shown in Figure 1.</figcaption></figure>

For $\mathcal{E}$ and $\mathcal{L}$, the equilibration is again somewhat slower for $\gamma = 0.5 \rightarrow 0.05$ than the reverse. We can further identify clear scalings $\mathcal{E} \sim D^{-1}$ and $\mathcal{L} \sim D^{-1/2}$. Finally, $\mathcal{L}$ is greater for $\gamma = 0.5 \rightarrow 0.05$ than the reverse. These results are all understandable in terms of the interpretation of $\mathcal{L}$ as the number of statistically-distinguishable states that the PDF evolves through: First, we recall from Figure 1 that $\gamma = 0.5 \rightarrow 0.05$ had consistently narrower PDFs than the reverse. Narrower PDFs means more distinguishable states, hence larger $\mathcal{L}$ for $\gamma = 0.5 \rightarrow 0.05$ than the reverse. The $\mathcal{L} \sim D^{-1/2}$ scaling has the same explanation; smaller $D$ yields narrower PDFs, hence larger $\mathcal{L}$.

The first panel in Figure 3 shows the previous quantities $\langle x \rangle$ and $\mathcal{L} \cdot D^{1/2}$, but now plotted against each other rather than separately against time. The behaviour is exactly as one might expect, with $\mathcal{L}$ growing more or less linearly with distance from the initial position. The right panel in Figure 3 shows the entropy (15), again as a function of $\langle x \rangle$ rather than time, to emphasize the cyclic nature of the two processes. The significance is indeed as claimed above, with more localized PDFs having smaller entropy values. Note how $\gamma = 0.5 \rightarrow 0.05$, which had the narrower PDFs, has lower entropy values than the reverse process. Note also how reducing $D$ by a factor of two, thereby making the PDFs narrower, causes the entire cyclic pattern to shift downward by an essentially constant amount.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) shows $\mathcal{L} \cdot D^{1/2}$ and (b) entropy, both as functions of $\langle x \rangle$. Solid lines denote $\gamma = 0.5 \to 0.05$, dashed lines the reverse. Numbers besides curves indicate $D = 0.01, 0.02, 0.04$. The arrows on the entropy plot are a reminder of the direction of inward/outward motion.</figcaption></figure>

Figure 4 shows how the standard deviation, skewness and kurtosis behave, again, as functions of $\langle x \rangle$ throughout the two processes. The heavy green lines also show the behaviour that would be expected if the time-dependent PDFs were always gamma distributions throughout their evolution. That is, if gamma distributions have $\langle x \rangle = \gamma$, $\sigma = \sqrt{\gamma D}$, skewness = $2\sqrt{D/\gamma}$ and kurtosis = $6D/\gamma$ (setting $\epsilon = 1$), then expressed as functions of $\langle x \rangle$, we would have $(\sigma/D^{1/2}) = \sqrt{\langle x \rangle}$, ($\text{skewness}/\sqrt{D} = 2/\sqrt{\langle x \rangle}$ and ($\text{kurtosis}/D = 6/\langle x \rangle$). As we can see, the $\gamma = 0.5 \to 0.05$ process follows these functional relationships reasonably well (especially for skewness and kurtosis), but for $\gamma = 0.05 \to 0.5$, all three quantities deviate substantially.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) $\sigma/D^{1/2}$, (b) (skewness/D$^{1/2}$) and (c) (kurtosis/D), as functions of $\langle x \rangle$. Solid lines denote $\gamma = 0.5 \to 0.05$, dashed lines the reverse. Numbers besides curves indicate $D = 0.01, 0.02, 0.04$. The heavy green curves are $\sqrt{\langle x \rangle}$, $2/\sqrt{\langle x \rangle}$ and $6/\langle x \rangle$, respectively, and indicate the behaviour expected for exact gamma distributions.</figcaption></figure>

Further evidence of significant deviations from gamma distribution behaviour is seen in Figure 5, showing the difference (13) directly. As expected from Figure 4, $\gamma = 0.05 \to 0.5$ has a much greater difference than $\gamma = 0.5 \to 0.05$. The second and third panels show how the PDFs compare with the equivalent gamma distributions having the same $\langle x \rangle$ and $\sigma$ values as the actual PDFs at that instant. The differences are clearly visible, especially for $\gamma = 0.05 \to 0.5$, but also for $\gamma = 0.5 \to 0.05$.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) shows the difference (13) between the actual PDF and the equivalent gamma distribution, as functions of $\langle x \rangle$. Solid lines denote $\gamma = 0.5 \to 0.05$, dashed lines the reverse, with arrows also indicating the direction of motion. The dots at $\langle x \rangle = 0.3$ for $\gamma = 0.05 \to 0.5$ and $\langle x \rangle = 0.1$ for $\gamma = 0.5 \to 0.05$, correspond to the other two panels: (b) compares the $\gamma = 0.05 \to 0.5$ PDF with its equivalent gamma distribution; (c) compares the $\gamma = 0.5 \to 0.05$ PDF with its equivalent gamma distribution. The actual PDFs in each case are solid (red), and the equivalent gamma distributions are dashed (blue). $D = 0.04$ for both sets.</figcaption></figure>

4.2. $D > \gamma$

We next consider the case $D > \gamma$, where we demonstrated above that stationary solutions cannot exist at all, because the time-dependent PDF can only ever get closer and closer to the gamma distribution singularity at the origin, but can never actually achieve it. To explore what does happen in this case then, we simply repeat the above procedure, except that there is now only an ‘inward’ process, and no reverse. That is, instead of $\gamma = 0.5 \to 0.05$, let us consider $\gamma = 0.5 \to 0$ (Throughout this section, we will also take $D = 10^{-3}$, to facilitate comparison with results in the next section. For $\gamma = 0$, of course, any $D$ is greater than $\gamma$.).

Figure 6 shows the resulting PDFs and how they approach ever closer to the origin, but never actually achieve the $x^{-1}$ blow-up that would be implied by Equation (4) for $a = \gamma/D = 0$. The peak amplitude simply increases indefinitely, as $t^{1/2}$. The widths correspondingly also decrease; the apparent increase is an illusion caused by the logarithmic scale for $x$. The dashed lines also show the equivalent gamma distributions, as before. Note how the difference becomes increasingly noticeable; in line with the fact that the equivalent gamma distribution is tending toward its singular behaviour as $\langle x \rangle$ decreases, but the actual PDFs must always have $p(0) = 0$.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** The initial condition is a gamma distribution with $\gamma = 0.5$, $\epsilon = 1$ and $D = 10^{-3}$; $\gamma$ is then switched to zero, and the solution is evolved according to Equation (3). Numbers besides curves indicate time, from the initial condition at $t = 0$ to the final time 1000. The dashed curves indicate the equivalent gamma distributions having the same $\langle x \rangle$ and $\sigma$.</figcaption></figure>
Figure 7 is the equivalent of Figure 2 and directly compares $\gamma = 0.5 \to 0$ here with the previous $\gamma = 0.5 \to 0.05$. We see that $\langle x \rangle$ starts out very similarly, but instead of equilibrating to 0.05, it now tends to zero as $t^{-1}$. $\mathcal{E}$ again starts out similarly, but ultimately tends to zero much more slowly, as $t^{-3}$ instead of exponentially. This $t^{-3}$ scaling for $\mathcal{E}$ has an interesting consequence for $\mathcal{L}$, namely that $\mathcal{L}$ does saturate to a finite value $\mathcal{L}_{\infty}$ (since $\int t^{-3/2} dt$ remains bounded for $t \to \infty$) even though the PDF itself never settles to a stationary state.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** As in Figure 2, (a) shows $\langle x \rangle$; (b) shows $\mathcal{E} \cdot D$; and (c) $\mathcal{L} \cdot D^{1/2}$. Solid lines denote $\gamma = 0.5 \to 0$ for $D = 10^{-3}$, dashed lines the previous $\gamma = 0.5 \to 0.05$ for $D = 0.01$. Note how the scalings of $\mathcal{E}$ and $\mathcal{L}$ with $D$ are still preserved even when $D$ is changed by a factor of 10.</figcaption></figure>

Figure 8 shows entropy, $\sigma$, skewness and kurtosis, so some of the results as in Figures 3 and 4. Entropy and $\sigma$ are again both good measures of how narrow the PDF is, becoming ever smaller as the peak moves toward the origin. Skewness and kurtosis seem to follow the expected gamma distribution relationship extremely well, even though we saw before in Figure 6 that the PDFs are actually different from gamma distributions. As $\langle x \rangle \to 0$, both skewness and kurtosis thus become indefinitely large.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** (a) Entropy, (b) $\sigma/D^{1/2}$, (c) (skewness / $D^{1/2}$) and (d) (kurtosis / $D$), as functions of $\langle x \rangle$, for the $\gamma = 0.5 \to 0$ calculation from Figure 6. The heavy green curves in the last three panels are $\sqrt{\langle x \rangle}$, $2/\sqrt{\langle x \rangle}$ and $6/\langle x \rangle$, respectively, and indicate the behaviour expected for exact gamma distributions.</figcaption></figure>

### 4.3. $Q \neq 0$

Finally, we turn to the Fokker-Planck Equation (11) with additive noise included and use it to explore the two questions that could not be addressed otherwise. First, how does a process like $\gamma = 0.5 \to 0$ then equilibrate to a stationary solution? Second, what does the reverse process $\gamma = 0 \to 0.5$ look like?

We will keep $D = 10^{-3}$ and $Q = 10^{-5}$ fixed throughout this section. Since the effective diffusion coefficients in (11) are $Dx^2$ and $Q$ (recall also the denominator of Equation (12)), this means that $Q$ is
dominant only within $x \le 0.1$; any stationary solutions with peaks much beyond that are effectively pure gamma distributions.

Figure 9 shows the same type of inward/outward process as before in Figure 1, only now switching $\gamma$ between 0.5 and 0.1. Comparing with Figure 1, we see that the dynamics are very similar, just with all the peaks considerably narrower, which is to be expected if $D = 10^{-3}$ rather than 0.02. The only other point to note is how the final peak in the left panel is lower than the previous peak at $\langle x \rangle = 0.2$, which is different from Figure 1, where $\gamma = 0.5 \to 0.05$ had peaks monotonically increasing throughout the entire evolution. The reason for the final peak here decreasing slightly is precisely the influence of $Q$ in this region; if this peak is now seeing just as much diffusion from $Q$ as from $D$, it is not surprising that it spreads out somewhat more and is correspondingly somewhat lower than a pure gamma distribution would be.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** (a) shows the result of switching $\gamma = 0.5 \to 0.1$, (b) $\gamma = 0.1 \to 0.5$, both at fixed $\epsilon = 1$, $D = 10^{-3}$ and $Q = 10^{-5}$. The initial (red) and final (blue) gamma distributions are shown as heavy lines. The three intermediate lines are when the time-dependent solutions have $\langle x \rangle = 0.2, 0.3, 0.4$. $L_\infty = 25$ on the left and 16 on the right.</figcaption></figure>

Figure 10 shows the fundamentally new case, namely switching $\gamma$ between 0.5 and zero. The inward process $\gamma = 0.5 \to 0$ is again very similar to either Figure 1 or 9. The only difference from Figure 6 is that the process does actually equilibrate to a stationary solution now, as given by Equation (12). The reverse process $\gamma = 0 \to 0.5$ is rather different though. The initial central peak now broadens far more than previously seen in Figures 1 and 9.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** (a) the result of switching $\gamma = 0.5 \to 0$, (b) $\gamma = 0 \to 0.5$, both at fixed $\epsilon = 1$, $D = 10^{-3}$ and $Q = 10^{-5}$. The initial (red) and final (blue) gamma distributions are shown as heavy lines. The four intermediate lines are when the time-dependent solutions have $\langle x \rangle = 0.1, 0.2, 0.3, 0.4$. $L_\infty = 35$ on the left and 9.5 on the right.</figcaption></figure>
One interesting consequence of this extreme broadening for $\gamma = 0 \rightarrow 0.5$ is on the total information length $\mathcal{L}_{\infty}$. In Figure 9, these values are 25 and 16, respectively, whereas in Figure 10, they are 35 and 9.5. That is, in both cases, decreasing $\gamma$ yields larger $\mathcal{L}_{\infty}$ values than increasing $\gamma$ does, consistent with the peaks being narrower, hence passing through more statistically-distinguishable states. Next, comparing $\mathcal{L}_{\infty} = 25$ for $\gamma = 0.5 \rightarrow 0.1$ versus 35 for $\gamma = 0.5 \rightarrow 0$, this is exactly as one might expect: having the peak travel somewhat further yields extra information length. However, comparing $\mathcal{L}_{\infty} = 16$ for $\gamma = 0.1 \rightarrow 0.5$ versus 9.5 for $\gamma = 0 \rightarrow 0.5$ is puzzling then. The peak has further to travel, but accomplishes it with less information length. The reason is precisely this extreme broadening, which substantially reduces the number of distinguishable states along the way. See also [40,41], where the same effect was studied for Gaussian PDFs and values of $D$ as small as $10^{-7}$, leading to fundamentally different scalings of $\mathcal{L}_{\infty}$ with $D$ for inward and outward processes.

Returning to the central question of this paper, namely how close the time-dependent PDFs are to gamma distributions, the results for Figure 9 are similar to the previous ones. In particular, we recall that before in Figure 5, we had the difference scaling as $D^{1/2}$, so a smaller $D$ here means a smaller difference. These results are approaching the small $D$ regime where gamma distributions become very close to Gaussians anyway, which generally remain close to Gaussian as they move.

However, for the $\gamma = 0 \rightarrow 0.5$ process in Figure 10, the intermediate stages do not look much like gamma distributions (the final equilibrium is indistinguishable from a gamma distribution though, consistent with $Q$ being completely negligible for these values of $x$). For the intermediate stages, these were found to be so different from gamma distributions that attempting to fit a gamma distribution having the same $\langle x \rangle$ and $\sigma$ made little sense; this extreme broadening and long tail trailing behind the peak meant that both $\langle x \rangle$ and $\sigma$ were too different from the normal expectation that they should be measures of ‘peak’ and ‘width’.

Instead, we simply asked the question, which values of $a$ and $b$ would minimize the quantity $\int |p - p_{bf}| dx$, where $p$ is the time-dependent PDF to be fitted and $p_{bf}$ is the best-fit gamma distribution. Unlike our previous difference formula, this does not yield simple analytic formulas for the $a$ and $b$ to choose, but is numerically still straightforward to implement. Figure 11 shows the results, for two of the intermediate stages in the $\gamma = 0 \rightarrow 0.5$ process. We can see that the fit is rather poor, indicating that these PDFs are significantly different from gamma distributions.

This misfit is also not caused by the inclusion of $Q$; if this or any similar central peak is evolved for either small or zero $Q$ in the Fokker-Planck equation, the result is always similar to here. As explained also in [40,41], the dynamics of how central peaks move away from the origin is simply different from how peaks already away from the origin move, regardless of whether the final states are Gaussians as in [40,41] or gamma distributions as here.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** The $\gamma = 0 \rightarrow 0.5$ process as in Figure 10, but now shown in more detail. The dashed (magenta) curves are the gamma distributions that best fit the two thicker curves at intermediate times. Note how even a “best-fit” is a rather poor approximation to the actual PDFs.</figcaption></figure>
## 5. Conclusions

Gamma distributions are among the most popular choices for modelling a broad range of experimentally-determined PDFs. It is often assumed that time-dependent PDFs can then simply be modelled as gamma distributions with time-varying parameters a and b. In this work, we have demonstrated that one should be cautious with such an approach. By numerically solving the full time-dependent Fokker-Planck equation, we found that there are three sets of circumstances where the PDFs can differ significantly from gamma distributions:

* If $D < \gamma$, so that stationary solutions exist, but D is also sufficiently close to $\gamma$ that a gamma distribution differs significantly from a Gaussian, then the time-dependent PDFs will also differ significantly from gamma distributions.

* If $D > \gamma$, stationary gamma distributions do not exist at all. Instead, peaks move ever closer to the origin and in the process increasingly differ from gamma distributions.

* If the initial condition is a peak right on the origin—either as a result of adding additive noise to produce stationary solutions even for $D > \gamma$, or simply as an arbitrary initial condition—then any evolution away from the origin will differ significantly from gamma distributions. Unlike the previous two items, which become more pronounced for larger D, this effect is most clearly visible for smaller D, where the mismatch between the naturally narrower peaks and the extreme broadening seen in Figure 11 becomes increasingly significant.

In summary, our results show that a simple Langevin equation model mimics the strong fluctuations of far-from-equilibrium systems. This model has gamma distributions as steady-state solutions, but the time-dependent solutions can deviate considerably from this law. This makes tasks such as Bayesian and frequentist inference of the model from data more complicated. On the other hand, the model shows complex asymptotic dynamics with situations when a steady state is reached or not, different from the one of the deterministic logistic model that invariably evolves to the maximum capacity. The studied model is general enough and can be applied to many practical situations in biology, economics, finance and physics. Future work will apply some of these ideas to fitting actual data.

**Author Contributions:** Eun-jin Kim and Rainer Hollerbach conceived the basic mathematical ideas; Ovidiu Radulescu provided background in biology and other areas; Eun-jin Kim conducted the analytical derivations; Lucille-Marie Tenkès and Rainer Hollerbach conducted the numerical calculations; all authors were involved in writing the paper. All the authors have read and approved the final manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A. Derivation of the Fokker-Planck Equations

In order to derive the Fokker-Planck Equation (3) from the Langevin Equation (1), it is useful to introduce a generating function Z:

$$Z = e^{i\lambda x(t)}. \qquad (A1)$$

Then, by definition of 'average', the average of Z is related to the PDF, $p(x, t)$, as:

$$\langle Z \rangle = \int dx Z p(x, t) = \int dx e^{i\lambda x(t)} p(x, t). \qquad (A2)$$

Thus, we see that $\langle Z \rangle$ is the Fourier transform of $p(x, t)$. The inverse Fourier transform of $\langle Z \rangle$ then gives $p(x, t)$:

$$p(x, t) = \frac{1}{2\pi} \int d\lambda e^{-i\lambda x} \langle Z \rangle. \qquad (A3)$$

We note that Equation (A3) can be written as:

$$p(x, t) = \left\langle \frac{1}{2\pi} \int d\lambda e^{i\lambda(x-x(t))} \right\rangle = \langle \delta(x-x(t)) \rangle, \qquad (A4)$$
which is another form of $p(x, t)$. To obtain the equation for $p(x, t)$, we first derive the equation for $\langle x \rangle$ and then take the inverse Fourier transform as summarised in the following.

We differentiate $Z$ with respect to time $t$ and use Equation (1) to obtain:

$$ \partial_t Z = i\lambda\partial_t xZ = i\lambda(\gamma x - \epsilon x^2 + \xi(t)x)Z = \lambda[\gamma\partial_\lambda + i\epsilon\partial_{\lambda\lambda} + \xi\partial_\lambda]Z, \quad (A5) $$

where $xZ = -i\partial_\lambda Z$ was used. The formal solution to Equation (A5) is:

$$ Z(t) = \lambda \int dt_1 [\gamma \partial_\lambda + i\epsilon \partial_{\lambda\lambda} + \xi(t_1) \partial_\lambda] Z(t_1). \quad (A6) $$

The average of Equation (A5) gives:

$$ \partial_t \langle Z \rangle = \lambda(\gamma\partial_\lambda + i\epsilon\partial_{\lambda\lambda})\langle Z \rangle + \lambda\langle\xi(t)\partial_\lambda Z(t)\rangle. \quad (A7) $$

To find $\langle \xi(t) \partial_\lambda Z(t) \rangle$, we use Equation (A6) iteratively as follows:

$$ \begin{aligned} \langle \xi(t) \partial_\lambda Z \rangle &= \left\langle \xi(t) \partial_\lambda \left[ \lambda \int dt_1 [\gamma \partial_\lambda + i\epsilon \partial_{\lambda\lambda} + \xi(t_1) \partial_\lambda] Z(t_1) \right] \right\rangle \\ &= \langle \xi(t) \rangle \partial_\lambda \left[ \lambda \int dt_1 [\gamma \partial_\lambda + i\epsilon \partial_{\lambda\lambda}] \langle Z(t_1) \rangle \right] + \partial_\lambda \left[ \lambda \int dt_1 \langle \xi(t) \xi(t_1) \rangle \partial_\lambda \langle Z(t_1) \rangle \right] \\ &= \partial_\lambda \left[ \lambda [D \partial_\lambda \langle Z(t) \rangle] \right]. \end{aligned} \quad (A8) $$

Here, we used the independence of $\xi(t)$ and $Z(t_1)$ for $t_1 < t$, $\langle \xi(t) Z(t_1) \rangle = \langle \xi(t) \rangle \langle Z(t_1) \rangle = 0$, together with Equation (2), $\int_0^t dt_1 \delta(t - t_1) = 1/2$, and $\langle \xi \rangle = 0$. By substituting Equation (A8) into Equation (A7), we obtain:

$$ \partial_t \langle Z \rangle = \lambda(\gamma\partial_\lambda + i\epsilon\partial_{\lambda\lambda})\langle Z \rangle + \lambda\partial_\lambda \left[ \lambda [D\partial_\lambda \langle Z(t) \rangle] \right]. \quad (A9) $$

The inverse Fourier transform of Equation (A9) then gives us:

$$ \frac{\partial}{\partial t} p(x, t) = -\frac{\partial}{\partial x} \left[ (\gamma x - \epsilon x^2) p(x, t) \right] + D \frac{\partial}{\partial x} \left[ x \frac{\partial}{\partial x} \left[ x p(x, t) \right] \right] \quad (A10) $$

which is Equation (3). Specifically, the inverse Fourier transforms of the first and last terms in Equation (A9) are shown explicitly in the following:

$$ \frac{1}{2\pi} \int d\lambda e^{-i\lambda x} \langle \partial_t Z \rangle = \partial_t \left[ \frac{1}{2\pi} \int d\lambda e^{-i\lambda x} \langle Z \rangle \right] = \frac{\partial}{\partial t} p(x, t), \quad (A11) $$

$$ \frac{D}{2\pi} \int d\lambda e^{-i\lambda x} \lambda \partial_\lambda [\lambda \partial_\lambda \langle Z \rangle] = D \frac{\partial}{\partial x} \left[ x \frac{\partial}{\partial x} \left[ x p(x, t) \right] \right], \quad (A12) $$

where integration by parts was used twice in obtaining Equation (A12). The additional $Q\partial_{xx}p$ term in the Fokker-Planck Equation (11) can be derived in the same way.

## Appendix B. Time-Dependent Analytical Solutions of Equation (3)

We begin by making the change of variables $y = 1/x$ in Equation (1) to obtain:

$$ \frac{dy}{dt} = -(\gamma + \xi)y + \epsilon. \quad (A13) $$

By using the Stratonovich calculus [2,3,38], the solution to Equation (A13) is found as:

$$ y(t) = y_0 e^{-(\gamma t + B(t))} + \epsilon e^{-(\gamma t + B(t))} \int_0^t dt_1 e^{(\gamma t_1 + B(t_1))}, \quad (A14) $$
where $y_0 = y(t = 0)$ and $B(t) = \int_0^t dt_1 \xi(t_1)$ is the Brownian motion. Therefore,

$$x(t) = \frac{x_0 e^{\gamma t + B(t)}}{1 + \epsilon x_0 \int_0^t dt_1 e^{(\gamma t_1 + B(t_1))}}, \quad (A15)$$

where $x_0 = x(t = 0)$. In Equation (A15), $e^{B(t)}$ is the geometric Brownian motion while $e^{-\gamma t - B(t)}$ is the geometric Brownian motion with a drift (e.g., [2]). The time integral of the latter is used in understanding stochastic processes in financial mathematics and many other areas [49,50]. In particular, in the long time limit, its PDF can be shown to be a gamma distribution. However, this PDF of $x$ is not particularly useful as it involves complicated summations and integrals that cannot be evaluated in closed form [49,50].

## References

1. Risken, H. *The Fokker-Planck Equation: Methods of Solution and Applications*; Springer: Berlin, Germany, 1996.

2. Klebaner, F. *Introduction to Stochastic Calculus with Applications*; Imperial College Press: London, Britain, 2012.

3. Gardiner, C. *Stochastic Methods*, 4th ed.; Springer: Berlin, Germany, 2008.

4. Saw, E.-W.; Kuzzay, D.; Faranda, D.; Guittonneau, A.; Daviaud, F.; Wiertel-Gasquet, C.; Padilla, V.; Dubrulle, B. Experimental characterization of extreme events of inertial dissipation in a turbulent swirling flow. *Nat. Commun.* **2016**, *7*, 12466, doi:10.1038/ncomms12466.

5. Kim, E.; Diamond, P.H. On intermittency in drift wave turbulence: Structure of the probability distribution function. *Phys. Rev. Lett.* **2002**, *88*, 225002.

6. Kim, E.; Diamond, P.H. Zonal flows and transient dynamics of the L-H transition. *Phys. Rev. Lett.* **2003**, *90*, 185006.

7. Kim, E. Consistent theory of turbulent transport in two dimensional magnetohydrodynamics. *Phys. Rev. Lett.* **2006**, *96*, 084504.

8. Kim, E.; Anderson, J. Structure-based statistical theory of intermittency. *Phys. Plasmas* **2008**, *15*, 114506.

9. Newton, A.P.L.; Kim, E.; Liu, H.-L. On the self-organizing process of large scale shear flows. *Phys. Plasmas* **2013**, *20*, 092306.

10. Srinivasan, K.; Young, W.R. Zonostrophic instability. *J. Atmos. Sci.* **2012**, *69*, 1633–1656.

11. Sayanagi, K.M.; Showman, A.P.; Dowling, T.E. The emergence of multiple robust zonal jets from freely evolving, three-dimensional stratified geostrophic turbulence with applications to Jupiter. *J. Atmos. Sci.* **2008**, *65*, 3947–3962.

12. Tsuchiya, M.; Giuliani, A.; Hashimoto, M.; Erenpreisa J.; Yoshikawa, K. Emergent self-organized criticality in gene expression dynamics: Temporal development of global phase transition revealed in a cancer cell line. *PLoS ONE* **2015**, *10*, e0128565.

13. Tang, C.; Bak, P. Mean field theory of self-organized critical phenomena. *J. Stat. Phys.* **1988**, *51*, 797–802.

14. Jensen, H.J. *Self-Organized Criticality: Emergent Complex Behavior in Physical and Biological Systems*; Cambridge University Press: Cambridge, UK, 1998.

15. Pruessner, G. *Self-Organised Criticality*; Cambridge University Press: Cambridge, UK, 2012.

16. Longo, G.; Montévil, M. From physics to biology by extending criticality and symmetry breaking. *Prog. Biophys. Mol. Biol.* **2011**, *106*, 340–347.

17. Flynn, S.W.; Zhao, H.C.; Green, J.R. Measuring disorder in irreversible decay processes. *J. Chem. Phys.* **2014**, *141*, 104107.

18. Nichols, J.W.; Flynn, S.W.; Green, J.R. Order and disorder in irreversible decay processes. *J. Chem. Phys.* **2015**, *142*, 064113.

19. Ferguson, M.L.; Le Coq, D.; Jules, M.; Aymerich, S.; Radulescu, O.; Declerck, N.; Royer, C.A. Reconciling molecular regulatory mechanisms with noise patterns of bacterial metabolic promoters in induced and repressed states. *Proc. Natl. Acad. Sci. USA* **2012**, *109*, 155–160.

20. Shahrezaei, V.; Swain, P.S. Analytical distributions for stochastic gene expression. *Proc. Natl. Acad. Sci. USA* **2008**, *105*, 17256.

21. Thomas, R.; Torre, L.; Chang, X.; Mehrotra, S. Validation and characterization of DNA microarray gene expression data distribution and associated moments. *BMC Bioinform.* **2010**, *11*, 576.

22. Iyer-Biswas, S.; Hayot, F.; Jayaprakash, C. Stochasticity of gene products from transcriptional pulsing. *Phys. Rev. E* **2009**, *79*, 031911.
23. Elgart, V.; Jia, T.; Fenley, A.T.; Kulkarni, R. Connecting protein and mRNA burst distributions for stochastic models of gene expression. *Phys. Biol.* **2011**, *8*, 046001.

24. Delbrück, M. The burst size distribution in the growth of bacterial viruses (bacteriophages). *J. Bacteriol.* **1945**, *50*, 131.

25. Delbrück, M. Statistical fluctuations in autocatalytic reactions. *J. Chem. Phys.* **1940**, *8*, 120–124.

26. Kim, E.; Liu, H.; Anderson, J. Probability distribution function for self-organization of shear flows. *Phys. Plasmas* **2009**, *16*, 052304.

27. Kim, E.; Hollerbach, R. Time-dependent probability density function in cubic stochastic processes. *Phys. Rev. E* **2016**, *94*, 052118.

28. De Jong, I.G.; Haccou, P.; Kuipers, O.P. Bet hedging or not? A guide to proper classification of microbial survival strategies. *Bioessays* **2011**, *33*, 215–223.

29. Balaban, N.; Merrin, J.; Chait, R.; Kowalik, L.; Leibler, S. Bacterial persistence as a phenotypic switch. *Science* **2004**, *305*, 1622–1625.

30. Glansdorff, P.; Prigogine, I. Thermodynamic theory of structure, stability and fluctuations. *Am. J. Phys.* **1973**, *41*, 147–148.

31. Suzuki, M. Microscopic theory of formation of macroscopic order. *Phys. Lett. A* **1980**, *75*, 331–332.

32. Suzuki, M. The variational theory and rate equation method with applications to relaxation near the instability point. *Phys. A Stat. Mech. Appl.* **1981**, *105*, 631–641.

33. Langer, J.S.; Baron, M.; Miller, H.D. New computational method in the theory of spinodal decomposition. *Phys. Rev. A* **1975**, *11*, 1417.

34. Saito, Y. Relaxation in a bistable system. *J. Phys. Soc. Jpn.* **1976**, *61*, 388–393.

35. Hasegawa, H. Variational approach in studies with Fokker-Planck equations. *Prog. Theor. Phys.* **1977**, *58*, 128–146.

36. Dennis, B.; Costantino, R.F. Analysis of steady-state populations with the gamma abundance model: Application to Tribolium. *Ecology* **1988**, *69*, 1200–1213.

37. Liao, H.-Y.; Ai, B.-Q.; Hu, L. Effects of multiplicative colored noise on bacteria growth. *Braz. J. Phys.* **2007**, *37*, 1125–1128.

38. Wong, E.; Zakai, M. On the convergence of ordinary integrals to stochastic integrals. *Ann. Math. Stat.* **1960**, *36*, 1560–1564.

39. Bagui, S.C.; Mehra, K.L. Convergence of binomial, Poisson, negative-binomial, and gamma to normal distribution: Moment generating functions technique. *Am. J. Math. Stat.* **2016**, *6*, 115–121.

40. Kim, E.; Hollerbach, R. Geometric structure and information change in phase transitions. *Phys. Rev. E* **2017**, *95*, 062107.

41. Hollerbach, R.; Kim, E. Information geometry of non-equilibrium processes in a bistable system with a cubic damping. *Entropy* **2017**, *19*, 268.

42. Tenkès, L.-M.; Hollerbach, R.; Kim, E. Time-dependent probability density functions and information geometry in stochastic logistic and Gompertz models. *arXiv* **2017**, arXiv:1708.02789.

43. Frieden, B.R. *Physics from Fisher Information*; Cambridge University Press: Cambridge, Britain, 2000.

44. Wootters, W.K. Statistical distance and Hilbert space. *Phys. Rev. D* **1981**, *23*, 357.

45. Nicholson, S.B.; Kim, E. Investigation of the statistical distance to reach stationary distributions. *Phys. Lett. A* **2015**, *379*, 83–88.

46. Nicholson, S.B.; Kim, E. Structures in sound: Analysis of classical music using the information length. *Entropy* **2016**, *18*, 258, doi: 10.3390/e18070258.

47. Heseltine, J.; Kim, E. Novel mapping in a non-equilibrium stochastic process. *J. Phys. A* **2016**, *49*, 175002.

48. Kim, E.; Lee, U.; Heseltine, J.; Hollerbach, R. Geometric structure and geodesic motion in a solvable model of non-equilibrium stochastic process. *Phys. Rev. E* **2016**, *93*, 062127.

49. Bertoin, J.; Yor, M. Exponential functionals of Lévy processes. *Prob. Surv.* **2005**, *2*, 191–212.

50. Matsumoto, H.; Yor, M. Exponential functionals of Brownian motion, I: Probability laws at fixed time. *Prob. Surv.* **2005**, *2*, 312–347.

<figure><img src="image_13.png" /></figure>