The Fokker–Planck equation of the superstatistical fractional Brownian motion with application to passive tracers inside cytoplasm

By collecting from literature data experimental evidence of anomalous diffusion of passive tracers inside cytoplasm, and in particular of subdiffusion of mRNA molecules inside live Escherichia coli cells, we obtain the probability density function of molecules’ displacement and we derive the corresponding Fokker–Planck equation. Molecules’ distribution emerges to be related to the Krätzel function and its Fokker–Planck equation to be a fractional diffusion equation in the Erdélyi–Kober sense. The irreducibility of the derived Fokker–Planck equation to those of other literature models is also discussed.


Introduction
The experimental evidence of anomalous diffusion in living systems has been definitively established [1,2,3] and, in particular, we remind the measurements of the motion of mRNA molecules inside live E. coli cells by Golding & Cox [4] that are now a milestone in the field.Here we derive the Fokker-Planck equation for the probability density function (PDF) of molecules' displacement in agreement with the main findings from such dataset and similars.
Unfortunately, the PDF of passive tracers in cytoplasm is not available yet from data because of technical issues that limit the number of particle's trajectories and this affects, in particular, the reliability of the tails of the distribution which are the footprint of deviation from Gaussianity and standard diffusion.Thus, in these circumstances with lack of measurements, we do not analyse experimental data but we collect experimental results [4,5,6,7,8,9,10,11] in order to derive first the molecules' PDF and then the corresponding Fokker-Planck equation.If the solution, namely the PDF, is already known then the governing equation can be considered useless but, in the following, we discuss how indeed knowing the PDF is not all in the game and the knowledge of the corresponding equation provides indeed futher valuable information.
The superstatistical fBm is indeed a randomly scaled Gaussian process that was studied for fractional anomalous diffusion originally within the framework of the generalised gray Brownian motion (ggBm) [29,30,31,32,28], that recently has been extended to investigate the relation with generalised time-fractional diffusion equations [33,34].Moreover, randomly scaled Gaussian processes resembling the ggBm have been formulated for modelling and understanding anomalous diffusion also within an underdamped approach [35,36,37].
The qualitative success of the superstatistical fBm for modelling anomalous diffusion in living systems was already discussed within the framework of the ggBm in relation to ergodicity breaking and fractional diffusion [28].The PDF of molecules is therefore a mixture of Gaussian distributions with random variances and, within a Bayesian approach, such population can be understood as the likelihood modulated by the prior distribution of a parameter.The formal randomization of this parameter is equivalent to the computation of the marginal likelihood, which corresponds indeed to the PDF of i.i.d.random variables [38].This view allows for highlighting and clarifing the role of the central limit theorem in the dynamics of an heterogeneous ensemble of Brownian particles as in the superstatistical fBm here considered [38].Any value of the variance of the Gaussian PDFs is idiosyncratic and the random diffusion coefficients represent the heterogeneity of the ensemble of molecules, which is sufficient to (weakly) break the ergodicity of the system [28].
To conclude, we report that the superstatical modelling of anomalous diffusion can be related also with the effect of probes' polidispersity that generates an apparent anomalous diffusion as emerged in the cytoplasm of human cells where the apparent anomaly exponent decreases with increasing polydispersity of the probes [39].These results can be applied also in intracellular studies of the mobility of nanoparticles, polymers, or oligomerizing proteins.
Here we found that the PDF of molecule displacement belongs to the family of the Krätzel special function [40] and its Fokker-Planck equation is a fractional diffusion equation in the Erdélyi-Kober sense [41,42] that cannot be reduced to existing models for anomalous diffusion.
The rest of the paper is organised as follows.In the next Section 2 we present the dataset by Golding & Cox [4] and the related findings from the corresponding data analysis.In Section 3 we establish the molecule PDF and in Section 4 we derive the governing Fokker-Planck equation.In Section 5 we discuss the results and in particular the valuable information behind the determination of the Fokker-Planck equation in spite of the known PDF.The last Section 6 is devoted to the conclusions and future perspectives.

The paradigmatic dataset by Golding & Cox
In this study, we are interested in the governing equation of the PDF of passive tracers diffusing in cytoplasm because of the information that it can give about the process, see §5, apart from the PDF itself that is its fundamental solution.In particular, for our aims, we mainly consider the results coming from the data acquired by Golding & Cox [4].The trajectories included in that dataset are 21 and so they are not enough for calculating ensemble averages.In this respect, we remind here that anomalous diffusion is indeed characterised also by a non-Gaussian distribution of particles and therefore also by the scaling of the tails of the distribution which are indeed affected by the size of the sample of the observations.
In their experiment, they considered the motion of mRNA molecules released from their template DNA and free to move in the cytoplasm, in particular, they tracked the random motion of individual fluorescently labeled mRNA molecules inside live E. coli cells.This dataset turned to be a benchmark dataset for studying anomalous diffusion in living systems and widely analysed in these years, see, e.g., [5,6,7,8,9,10,11].The main findings from the Golding & Cox experiment [4] can be summerized as follows.
Let X t be the molecule position at time t > 0, the time-averaged mean square displacement (TA-MSD) of molecules resulted to be subdiffusive [X t+∆ − X t ] 2 ∼ ∆ β [4], with β ∈ (0, 1).Moreover, TA-MSD emerged to be quite scattered [4,10,11] and this can be due to a population of diffusion coefficients that leads to weak ergodicity breaking [18,28].Furthermore, the application of the method of p-variation showed that the underlying motion of the molecules is more likely the fBm [5,6,10].Therefore, by adopting a superstatistics of fBm, namely by adopting the process X t = √ Λ B H t where Λ is a non-negative random variable and B H t is the fBm with Hurst exponent H ∈ (0, 1), it was possible to estimate the diffusion coefficients of the molecules' ensemble which resulted in a population related to the Weibull distribution [10], that actually is a special case of the generalised Gamma distribution.The ggBm is recovered when Λ is a totally skewed positive α-stable random variable with stable index α ∈ (0, 1) [29,28].A Weibull distribution of the diffusion coefficients has been confirmed [11] also when the fBm was replaced by the so-called fractional Lévy stable motion (FLSM) that is a Lévy-driven stochastic process that generalizes the fBm [7,11].
Trajectories from the Golding-Cox dataset are not enough for properly testing mixing and ergodicity by comparing ensemble and time averages.However, necessary conditions for mixing and ergodicity can be derived on the basis of a single trajectory [9].Trajectories from the Golding-Cox dataset satisfy the necessary conditions for mixing and ergodicity [9], which is consistent with both the fBm [10] and the proper FLSM model [8,9].Actually, if each trajectory is properly re-scaled then their scattering is largely reduced and the ergodicity breaking is no longer present [10].
Since the FLSM is based on Lévy stable distributions, diverging statistics as, for example, the MSD are replaced by the sample MSD that may exhibit either normal and anomalous diffusion.The FLSM results to perform comparable or even better than the fBm, at least for some of the trajectories, for what concerns some observable as the Hurst exponent, stability index and both sample MSD and sample p-variation [7].This may be due to the fact that the experimental trajectories may deviate from Gaussianity, contrary to the fBm approach, and then the more flexible FLSM allows for catching this feature [11].
Actually, further valuable experimental data on diffusion of passive tracers in cytoplasm have been published and so their analysis, too, beside those related with the Golding-Cox dataset [4].Here, we remind, for example, the studies of particle motion in crowded fluids, that is an analogue of the motion inside the cytoplasm of living cells, by using dextran dissolved in water [43,44,45] or purely viscous solution obtained by using sucrose into water [45] that led to the understanding of the corresponding diffusive motions in terms of the fBm [20,22], as well as the measurements of the dynamics of histonelike nucleoid-structuring proteins in live E. coli bacteria [46] where a power-law distribution of the diffusion coefficients of individual proteins has been observed in agreement with the Pearson Type VII distribution that led to the the development of a modelling approach based on the fBm that hierarchically takes into account the joint fluctuations of both the anomalous diffusion exponents and the diffusion constants [24], or diffusion of tracer particles in the cytoplasm of mammalian cells where the experimental observations are described by an intermittent fBm alternating between two states of different mobility [22] and similarly in the case of intracellular endosomes [23] suggesting that the underlying trajectories can be modelled by a fBm with a distributed Hurst exponent [25].

Superstatistical molecules' distribution
In the case under consideration, the PDF of molecules results to be given by the superstatistical integral where is the Guassian PDF of the fBm and f (λ) is, in short notation, the generalised Gamma distribution with λ, λ 0 , ν, ρ > 0. The generalised Gamma distribution reduces to the Weibull distribution when ν = ρ, i.e., f (λ; λ 0 , ρ, ρ) = W (λ; λ 0 , ρ), to the Gamma distribution when ρ = 1 and to the exponential distribution when ρ = ν = 1, see a comparison in Figure 1.The generalised Gamma distribution was studied in the framework of superstatistics both in the original formulation [47] and in the recent formulations for anomalous diffusion within the diffusing-diffusivity [48,49] and the ggBm approaches [48].
We observe that molecules' PDF (1) can be re-written as where and the variance of particle displacement is Plots of the molecules' PDF (1) are shown in Figure 2.Moreover, we consider now the Mellin transform pair [50] M r {ϕ; s} where L is a specific contour path that separates the poles which provide through the residue theorem a power-series with positive exponents from those poles that provide a power-series with negative exponents.Thus we obtain the following Mellin-Barnes integral representation For mathematical convenience, we introduce now the change of variable u = x 2 / 4λ 0 t 2H and, from the normalization constraint P (x; t) dx = P(u) du = 1, we have By using formula (9), the representation in terms of the H-Fox function [41,40] of the PDF (10) results to be and by considering the properties of Mellin-Barnes integrals [50], and so the properties of the H-Fox functions as well [51], we can get from (9) the following asymptotic expansions for u → 0 with d = min [−1/2, ν − 1], and for u → +∞ 4 Derivation of the Fokker-Planck equation The problem to derive the Fokker-Planck equation of ( 1) is similar to the problem to derive the relation between generalised diffusion equations and subordination schemes [52].Here, in order to derive the Fokker-Planck equation governing PDF (1), we first consider the following diffusion equation solved by the Gaussian PDF of the fBm (2) and, by multiplying both sides times f (λ) and integrating over λ, from formula (1) we have where the extended notation P (x; t, ν, ρ) for PDF (1) is used for highlighting the difference between the PDFs from both sides of the equation.
As a concluding remark, we compare the Fokker-Planck equation ( 18), or its special case (19), against other equations used in anomalous diffusion.From formula (17) the following noteworthy identities can be derived [41] where t D µ RL is the fractional derivative in the Riemann-Liouville sense [41], and we observe that formula (20) highlights the relation between Erdélyi-Kober fractional equations and Riemann-Liouville (or Caputo with proper initial conditions) fractional differential equations with time-varying coefficients [58].Thus, when ν = ρ = 1/(2H), equations ( 18) and (19) become that is different from existing fractional diffusion models with time-dependent diffusion coefficient [59,58,60,61] and, moreover, it cannot be reduced further up to the so-called time-fractional diffusion equation Equation ( 22) is the fractional Fokker-Planck equation, for example, of a continuous-time random walk (CTRW) with a fat-tailed distribution of waiting times, see, e.g., [62], but also for the intermediate asymptotic regime of a CTRW with two Markovian hopping-trap mechanisms [63].Such CTRW models have been used for modelling anomalous diffusion in living systems, see, e.g., subdiffusion of lipid granules in living fission yeast cells [17], and also, see, e.g., references [64,65,66], for explaining some statistical features that appear also in the Golding & Cox dataset [4].
The failure of the CTRW approach for modelling some common features of anomalous diffusion in living systems was already pointed-out [18,5,6] and an alternative and promising approach seemed to be the ggBm [28,10].As a matter of fact, the generalised Fokker-Planck equation for the ggBm emerged to be, as well, a fractional diffusion equation in the Erdélyi-Kober sense as follows [67]: but it is not related neither to (18) nor to (19) in anyone of their special cases.On the contrary, from (20) with ρ = 1/(2H), equation (23) reduces to (22), such that the governing equation of the CTRW is indeed a special case of the governing equation of the ggBm.The fundamental solution of ( 23) is provided by plugging in (1), in place of the generalised Gamma distribution (3), the following population of diffusion coefficients with 0 < 1/ρ < 1, where M β (z) is the M-Wright Mainardi function [68,69] and, as a matter of fact, the fundamental solution of (23) here denoted by P ggBm (x; t) is also an M-function [29,67], i.e., and then by setting ρ = 1/(2H) we have also the fundamental solution of ( 22) that is here denoted by P CTRW (x; t), i.e., The comparison between the generalised Gamma distribution of the diffusion coefficients f (λ) (3) and the M-function distribution f ggBm (λ) ( 24) is shown in Figure 3.Moreover, the comparison among the fundamental solutions (1), ( 25) and ( 26) of the the corresponding Fokker-Planck equations (18,19,21), ( 23) and ( 22), respectively, is shown in Figure 4 together with the Gaussian distribution.

Discussion
In spite of the fact that knowing the solution is the most desiderable observable, from the equation we can indeed obtain remarkable further information on the process.These information can be physical, mathematical and also useful for applications.
In fact, from the physical point-of-view we have already observed that the heterogeneity of the diffusion coefficients is the cause of the weak ergodicity breaking but from the Fokker-Planck equation we have also information about relation between physical quantities like, for example, the formula of the flux q(x, t), i.e., that, in the cosidered case, from the Fokker-Planck equation (18) it results [41] q(x, t) with n − 1 < 1/ρ ≤ 1, that is not proportional to the gradient as in standard diffusion and we have indeed that the heterogeneity is the cause of the emergence of a memory kernel.This memory kernel is determined by the specific distribution of the diffusion coefficients f (λ).Moreover, because of this memory kernel it follows that the resulting evolution equation is a fractional differential equation, and a special connection exists between the generalised Gamma distribution and the Erdélyi-Kober fractional operators, see [53,54,42].
Moreover, from the governing equation, we can see if there are forcing terms and their origin and, actually, we have from ( 18) that the heterogeneity of the diffusion coefficient causes a flux with memory but it does not cause the emerging of any apparent forcing: in fact the resulting governing equation does not include any terms like ∂F/∂x and preserves indeed the form of the free-particle diffusion equation with the feature of a timedependent effective diffusion coefficient.In particular, with respect to space, this heterogenity does not cause any effect and the operator in space remains the second derivative ∂ 2 /∂x 2 as in standard diffusion.
Furthermore, the derivation of Fokker-Planck equation ( 18) embodies as a matter of fact a derivation on physical ground of a fractional equation by avoiding the replacement of the operators.Actually, this derivation, together with the derivation of the governing equation of the ggBm [67], is an alternative derivation with respect to that on statistical grounds discussed in literature [53,54,42] of fractional differential equations in the Erdélyi-Kober sense: this physical significance of fractional equations let indeed to overpass questioned issues in fractional-dynamic generalizations and some related constraints that have to be checked [70].
From the mathematical point-of-view, we observe indeed that the process turns into a process governed by an integro-differential equation but it remains linear and parabolic, namely the emergence of a memory kernel is not the cause, for example, of the emergence of a second derivative in time, i.e, ∂ 2 P/∂t 2 , or a finite front velocity, by keeping the main characteristics of the standard diffusion.
To conclude this section about the importance to know the governing equation of a process rather than the solution in a special case, we report that, under the practical point of view and in the direction of applications, to know the equation allows for studying, at least numerically, the solution with particular initial and boundary conditions and also with real settings in complex and multidimensional domains, as well as, in presence of external forcings.

Conclusions and future perspectives
In this paper we focused on the signature of anomalous diffusion as measured in many biological systems and in particular as it is measured in the motion of passive tracers inside cytoplasm as the mRNA molecules inside live E. coli cells.From the wide literature, we paied special attention to the dataset by Golding & Cox [4] that has been largely analysed, see, e.g., [4,5,6,7,8,9,10,11].Among the many distinctive features, the most important for the present study is the determination of the single-molecule trajectory as a fBm, see, e.g., [5,19,6,20,21,22,23,24,25], and also the heterogeneity of the diffusion coefficient according to a Weibull-type distribution that emerges from the data analysis when the observed anomalous behaviour is reproduced through a fBm as the underlying stochastic process [10].
Following these experimental evidences, we slightly generalised the distribution of the diffusion coefficients by using a generalised Gamma distribution in analogy with other similar analysis [26,27,18] and we showed that the PDF of the mRNA molecule displacement is related with the Krätzel function and we derived the corresponding Fokker-Planck equation.The governing equation results to be a fractional diffusion equation in the Erdélyi-Kober sense with a time-dependent effective diffusion coefficient, which cannot be reduced to the governing equations of existing literature models as the CTRW or the ggBm, neither to other models governed by a fractional differential equation with a time-dependent diffusion coefficient.
This last fact is indeed a two-fold findings for motivating future investigations on this model setting.In fact, the derived equation ( 18) pushes the research towards the mathematical analysis of a novel family of equations as well as towards the development of solid and reliable numerical schemes for studying more realistic systems in multidimensional domains with real geometries and general initial and boundary condition, but also towards the development of statistical methods and tools for proper taking into account a distribution of diffusion coefficients that can lead to different modelling approaches as those based on an heterogeneous ensemble of particles, e.g., the over-and under-damped ggBm [29,28,35,48,10,36,71,37], or those based on diffusion in inhomogeneous random environments, e.g., the diffusing diffusivity approach [72,73,48,74,75].
Finally, we highlight that the present framework does not include, yet, two quite general and well established features of anomalous diffusion: the Brownian yet non-Gaussian regime [72,73,76,48,49,77,78] and the anomalous-to-normal transition [79,76,80,36] that are indeed part of the paradigm of anomalous diffusion [63] and observed both in the continuousspace setting of the diffusion diffusivity models [73,48] and in the discretespace setting of random walks [79] even with solely two Markovian hoppingtrap mechanisms [63].Thus, these last embody the future developments of research on the superstatistical fBm.