Far-From-Equilibrium Time Evolution between Two Gamma Distributions

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.

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,2,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,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,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 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,20,21,22,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,31,32,33,34,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 +\xi )x-ϵx^2,$$

(1)

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

$$\langle \xi (t)\xi (t^{\prime })\rangle =2D\delta (t-t^{\prime }).$$

(2)

In Equation (2), the angular brackets represent the average over $\xi$, $\langle \xi \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 $ϵ$ represents the efficiency in self-regulation by a negative feedback. $\gamma x-ϵx^2$ can be considered as the gradient of the potential V as $\gamma x-ϵx^2=-\frac{\partial V}{\partial x},$ where $V=-\frac{\gamma }{2}{x}^2+\frac{ϵ}{3}{x}^3$. Thus, V has its minimum value at $x=\frac{\gamma }{ϵ}$. When $\xi =0$, $x=\frac{\gamma }{ϵ}$ (the carrying capacity) is a stable equilibrium point since ${\partial }_{xx}V{|}_{x=\gamma /ϵ}=\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 $\xi$. 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-ϵx^2)p(x,t)\right]+D\frac{\partial }{\partial x}\left[x\frac{\partial }{\partial x}\left[xp(x,t)\right]\right]$$

(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},$$

(4)

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

$$\langle x\rangle =\frac{a}{b}=\frac{\gamma }{ϵ},\mathrm{Var}(x)={\sigma }^2=\langle {(x-\langle x\rangle )}^2\rangle =\frac{a}{b^2}=\frac{\gamma D}{{ϵ}^2},$$

(5)

where $\sigma =\sqrt{\mathrm{Var}(x)}$ is the standard deviation. We recognise $\langle x\rangle$ as the carrying capacity for a deterministic system with $\xi =0$. It is useful to note that $\langle x\rangle$ is given by the linear growth rate scaled by $ϵ$, while $\mathrm{Var}(x)$ is given by the product of the linear growth rate and the diffusion coefficient, each scaled by $ϵ$. 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 $ϵ$. In such a limit, a and b are large, making $\sqrt{\mathrm{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)},$$

(6)

where $f(x)=bx-(a-1)lnx$. 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:

$$\begin{array}{c}{x}_p=\frac{a-1}{b}\sim \frac{a}{b},f(x=x_p)\sim a\left(1-ln\frac{a}{b}\right),\hfill \end{array}$$

(7)

$$\begin{array}{c}f(x)\sim f(x_p)+\frac{1}{2}{(x-x_p)}^2{\partial }_{xx}f(x){|}_{x=x_p}=a\left(1-ln\frac{a}{b}\right)+\frac{b^2}{2a}{\left(x-\frac{a}{b}\right)}^2.\hfill \end{array}$$

(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],$$

(9)

which is a Gaussian PDF with mean value $\langle x\rangle$. Here, $\beta =1/\mathrm{Var}(x)$ is the inverse temperature, and the variance $\mathrm{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 ${\displaystyle \underset{x\to 0}{lim}}p={\displaystyle \underset{x\to 0}{lim}}\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 +\xi )x-ϵx^2+\eta ,$$

(10)

where $\xi$ and $\eta$ are uncorrelated and $\eta$ satisfies $\langle \eta (t)\eta (t^{\prime })\rangle =2Q\delta (t-t^{\prime })$. The new version of the Fokker–Planck Equation (3) then becomes:

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

(11)

which has stationary solutions given by:

$$lnp(x)=\int \frac{(\gamma -D)x-ϵx^2}{D{x}^2+Q}dx.$$

(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,41,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 pdx=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_{\mathrm{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:

$$\mathrm{Difference}=\int |p-p_{\mathrm{equiv}}|dx$$

(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:

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

(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 plnpdx.$$

(15)

In particular, we expect S to be small for localised PDFs and large for spread out ones (e.g., [40,41,42,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-lnb+ln(\Gamma (a))+(1-a)\psi (a),$$

(16)

where $\psi (a)=\frac{dln(\Gamma (x))}{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}{{\left[\tau (t)\right]}^2}=\int \frac{1}{p(x,t)}{\left[\frac{\partial p(x,t)}{\partial t}\right]}^{2}dx.$$

(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 $ϵ$ 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{d{t}_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}d{t}_1.$$

(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,41,42,43,44,45,46,47,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 $ϵ$ 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 $ϵ=1$ and $D=0.02$ (one of the three parameters $ϵ$, $\gamma$ and D can of course always be kept fixed by rescaling the entire equation, so throughout this entire section, we keep $ϵ=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].

For $\mathcal{E}$ and $\mathcal{L}$, the equilibration is again somewhat slower for $\gamma =0.5\to 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\to 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\to 0.05$ had consistently narrower PDFs than the reverse. Narrower PDFs means more distinguishable states, hence larger $\mathcal{L}$ for $\gamma =0.5\to 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}·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\to 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 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 $ϵ=1$), then expressed as functions of $\langle x\rangle$, we would have $(\sigma /D^{1/2})=\sqrt{\langle x\rangle }$, (skewness $/\sqrt{D})=2/\sqrt{\langle x\rangle }$ and (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.

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$.

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 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 8 shows entropy, $\sigma$, skewness and kurtosis, so some of the results as in Figure 3 and Figure 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.

4.3. $Q\ne 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 $D{x}^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 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 Figure 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 Figure 1 and Figure 9.

One interesting consequence of this extreme broadening for $\gamma =0\to 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\to 0.1$ versus 35 for $\gamma =0.5\to 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\to 0.5$ versus $9.5$ for $\gamma =0\to 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\to 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_{\mathrm{bf}}|dx$, where p is the time-dependent PDF to be fitted and $p_{\mathrm{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\to 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.

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.