Brownian non-Gaussian polymer diffusion and queing theory in the mean-field limit

We link the Brownian non-Gaussian diffusion of a polymer center of mass to a microscopic cause: the polymerization/depolymerization phenomenon occurring when the polymer is in contact with a monomer chemostat. The anomalous behavior is triggered by the polymer critical point, separating the dilute and the dense phase in the grand canonical ensemble. In the mean-field limit we establish contact with queuing theory and show that the kurtosis of the polymer center of mass diverges alike a response function when the system becomes critical, a result which holds for general polymer dynamics (Zimm, Rouse, reptation). Both the equilibrium and nonequilibrium behaviors are solved exactly as a reference study for novel stochastic modeling and experimental setup.


Introduction
Occurrence of Brownian yet non-Gaussian diffusion is increasingly reported in experiments analyzing thermally driven motion in complex biological contexts. Some examples are: (i) beads diffusing on lipid tubes [1], in networks [2,3], or in a matrix of micropillars [4]; (ii) tracers in colloidal, polymeric and active suspensions [5,6]; (iii) lipid molecules or proteins embedded in protein-crowded lipid membranes [7,8], in biological cells [9,10,11,12], in narrow corrugated channels with fluctuating cross-section [13], in anisotropic liquid crystals [14]; (iv) motion of individuals in heterogeneous populations such as nematodes [15]; (v) colloids in random force fields [16]. This interesting phenomenon motivates mesoscopic modeling invoking superposition of statistics [17,18,15,2], diffusing diffusivities [19,20,21,22,23,24,25,26], subordination concepts [20], continuous time random walk [27], and diffusion in disordered environments [28], but also calls for microscopic foundation, as the identification of underlying universal mechanisms could assist the interpretation of experimental evidence or inspire new protocols. In order to fill the latter gap, we recently proposed [29] that monomers aggregation and disaggregation in the polymerization/depolymerization process offers a natural foundation to such anomalous diffusion if one is monitoring the motion of the center of mass (CM), an idea which has been further explored in Ref. [30] through a many-body approach. Profiting of this intuition, here we put forward a most transparent universal phenomenon triggering Brownian yet non-Gaussian diffusion: a polymer in contact with a chemostatted monomer bath. By changing the monomer concentration in the bath, the polymer undergoes a shift from finite to infinite average growth, a contingency which separates the dilute from the dense polymer phase [31,32,33,34]. At the transition point, diverging size fluctuations trigger an initial non-Gaussian diffusion of the polymer CM which finally crosses over to an ordinary diffusion dynamics. We address this intriguing mechanism in the simplest possible contest in which chain polymerization [35] occurs as a birth-death process, and discuss it in the mean-field limit. By establishing contact with queuing theory, this enables us to obtain an explicit solution of the model in which the (short-time) kurtosis of the polymer CM, measuring the degree of non-Gaussianity, diverges as the polymerization process becomes critical. We emphasize that the solution is general enough to deal with the several polymer models known in the literature [32,36] (Rouse, Zimm, reptation), each giving rise to a different leptokurtic probability density functions (PDF). Moreover all basic, time-dependent statistical quantities can be explicitly derived. We finally note that, while the crossovertime to Gaussian diffusion is affected by critical slowing down, the reaction rate of the polymerization process still offers an independent parameter controlling how fast normal diffusion is restored.
The paper is organized as follows. In Section 2 we introduce the grand canonical partition function of chemostatted polymers at equilibrium and the mean-field master equations describing the stochastic polymerization/depolymerization process whose exact solution is detailed in Appendix B. In Section 3 we study the stochastic motion of the CM of chemostatted polymers and highlight its non-Gaussian behavior by computing exactly the time evolution of the kurtosis under both equilibrium and nonequilibrium conditions, together with the shape of the initial non-Gaussian PDF under equilibrium condition. These exact results are compared with Gillespie-Langevin simulations in Appendix D. A summary and discussion is provided in Section 4, while further mathematical details are deferred to the appendixes.

Critical polymers
Chemostatted polymers are conveniently described in the grand canonical ensemble where the monomer fugacity z governs the grand canonical partition function Z gc and the equilibrium distribution P N (n) associated to the event N = n for the fluctuating polymer size N . Close to criticality, z → z − c , the former behaves asymptotically as [31,32,33,34] where µ c is the (model-dependent) connective constant and z c = µ −1 c . The universal entropic exponent γ is specified by the space dimension d, by the underlying topology of the polymeric structure, and by the equilibrium phase: good, Θ-, or bad solvent (see, e.g., [37] and references therein). The critical point z = z c separates the dilute phase, characterized by a finite average size E[N ], from the dense one in which the average size diverges.
Interestingly, the equilibrium distribution can be related to a simple master equation describing the polymerization/depolymerization process occurring as monomers add and detach to the polymer in the grand canonical ensemble; in this way connection with queing theory is established. For convenience, if n min is the minimal polymer size in the chain polymerization process, we operate the change of variable ‡ n → n − n min which associates to P N (n) the support 0 ≤ n < ∞ without altering the asymptotic singular behavior close to criticality. Consider then the (forward) master equation ∂ t P N (n, t|n 0 ) = µ P N (n + 1, t|n 0 ) + λ(n − 1) P N (n − 1, t|n 0 ) −(µ + λ(n))P N (n, t|n 0 ) (n > 0) Here P N (n, t|n 0 ) is the probability for N = n at time t ≥ 0 given N = n 0 at t = 0, and λ(n), µ are the rates for association and dissociation, respectively. In ‡ Because of specificities of the mean-field limit we prefer this change of variable with respect to the one adopted in [37].
Eq. (3) it is assumed size-dependency for the association rate (as it typically relies on the local concentration of available monomers), whereas dissociation normally occurs independently of n. Defining the growth factor as g(n) ≡ λ(n)/µ, in the Appendix we show that stationarity is attained under detailed balance, g(n) = P N (n + 1)/P N (n): this identifies the polymerization process, given P N (n). Note that the rate µ remains a free parameter which may rescale Eq. (3), thus determining the time scale τ for the autocorrelation of N (t) (see below). Knowledge of γ is provided by a mapping to the magnetic O(n → 0) model [32]; in the present paper we focus on the mean-field limit, in which γ → 1. This provides the simplification λ(n) = λ , g(n) = g = z/z c ∀n .
In practical terms, the simplest situation we may associate to the mean-field limit is that of a linear polymer composed by N = n subunits A N , subject to the chemical reaction [38,35] where k + , k − are the rate constants for association and dissociation, respectively, and λ = k + c with c the (n-independent) local concentration of monomers A 1 . In this case, the grand canonical partition function and the equilibrium distribution recast into with n = −1/ ln g and U(n) the (discrete) unit step function. As g → 1 − , the exponential distribution widens tending to become uniform (over an infinite support). Consistently with what anticipated, we thus appreciate that g = 1 separates a g > 1 phase of infinite growth from a g < 1 phase of finite average size and variance respectively. Observe that for g 1 one may write σ 2 N = µ N /(1 − g) ∼ µ N since 1 − g is finite and g infinitesimal. On the other hand, as g → 1 − , the variance crosses over to the behavior σ 2 N = µ 2 N /g ∼ µ 2 N ; the latter being a signature of non-normal fluctuations arising at the critical point.
The mean-field master equation can be written as and we see that in this case the polymerization process corresponds in fact to the M/M/1 model (Markovian interarrival times/Markovian service times/1 server) in queuing theory [39]. The time-dependent solution of Eq. (7) is recapitulated in Appendix B, where it is also reported the time scale τ of the exponential decay of the auto-correlation coefficient, ρ(t) e −t/τ : Note that rescaling time by τ , renders the model independent of the reaction rate µ. The asymptotic behavior for small and large time of P N (n, t|n 0 ) is P N (n, t|n 0 ) ∼ t τ δ n,n 0 , P N (n, t|n 0 ) ∼ t τ P N (n), respectively.

Brownian non-Gaussian diffusion
A basic idea to provide diffusing diffusivity models with a microscopic footing is very simple. From polymer physics it is known that the CM position R CM of a macromolecule with N + n min subunits diffuses with a coefficient D(N ) = D 0 /(N + n min ) α , D 0 being specific to the subunit, and α to the chosen polymer model [32,36]. Notably, if the chain undergoes polymerization then N = N (t) becomes a stochastic process, and so does the diffusion coefficient D. In slow nucleation processes n min corresponds to the size of a polymer nucleus [40], and to keep contact with the exactly-solvable M/M/1 model the association and dissociation rates must be size-independent. By focusing on a linear polymer which can grow and deteriorate at both ends, in the following we consider n min =3. This means to assume the trimer as the minimal polymer conformation in the system. In this way, monomer can always attach and detach at the two extremities and λ, µ do not depend on the polymer size. Zimm model includes hydrodynamic interactions and via the Stokes-Einstein relation [36,41] the diffusion coefficient is proportional to the inverse of the polymer hydrodynamic radius R, D(N ) ∼ 1/R(N ) ∼ 1/N ν ; hence, α coincides with the mean-field metric exponent ν, α = ν = 1/2. Rouse dynamics is instead characterized by α = 1, and for reptation α = 2. On the Smoluchowski time scale §, R CM evolves according to where B(t) = N (0, t) is a Wiener process (Brownian motion) with infinitesimal increments dB(dt) = N (0, dt) -The notation N (µ, σ 2 ) indicates a Gaussian random variable with average µ and variance σ 2 . The same behavior is observed for any tagged monomer above the Rouse relaxation time [32,36]. It is important to stress that in Eq. (9) two sources of randomness are present. One is the standard thermal agitation imparted by the solvent and represented by B(t), the other is the polymerization process which affects the intensity of B(t) through the N -dependence of the diffusion coefficient. Technically, B(t) is referred to as the "subordinated process" and N (t) as the "subordinator process". Under ordinary conditions, the stationary distribution of § The Smoluchowski time scale is appropriate for the description of the polymer dynamics at any sizes N varying from single monomer to large colloids. Indeed, in the time needed to loose memory of inertial effects, the traveled distance with respect to the polymer radius a is given by In water at room temperature, this ratio varies from 10 −3 for single nucleotides or amino acids to 10 −5 for large colloids. N is strongly peaked around its mean value and the "diffusion of diffusivities effect" is then difficult to detect. The situation drastically changes in proximity of the critical point governing the divergence of the polymerization degree. For any given realization [n(t)] ≡ {n(t ) ∈ N | 0 ≤ t ≤ t} of the stochastic process N (t), the PDF of the CM location satisfies an ordinary diffusion equation: where r 0 is the initial CM position. Eq. (10) emphasizes the fact that, to determine p R CM (r, t|[n(t)]; r 0 ) at a given time t > 0, the knowledge of the whole history of the subordinator process is required. Only in this way at each update dt the correct diffusion coefficient can be provided to propagate the initial condition p R CM (r, 0|n 0 ; r 0 ) = δ(r−r 0 ) up to time t. By exploiting the scaling property c N (µ, σ 2 ) = N (c µ, c 2 σ 2 ) for c ∈ R, the diffusing path is conveniently reparametrized in terms of a coordinate which converts Eq. (9) into a standard overdamped Langevin equation with unit diffusion coefficient, where s ≥ 0 is a path variable corresponding to the realization of the stochastic process With respect to the "random path" s, Eq. (10) formally transforms into an ordinary diffusion equation [43], with Green function solution For a polymer of size n 0 + 3 starting with certainty at the origin, p R CM (r, 0|n 0 ; 0) = δ n,n 0 δ(r), the PDF of finding its CM at position r at time t is thus given by the subordination [44,45] formula where p S (s, t|n 0 ) is the probability distribution of the process S(t). Eq. (15)  The short time expansion of P N (n, t|n 0 ) reported in the Appendix exhibits nonlinear corrections (alternating series) to the linear increase of the mean squared displacement. The Brownian character becomes distinctive of large time, t τ , when P N (n, t|n 0 ) P N (n) and Eq. (16) simplifies to: with D av ≡ ∞ n=0 D 0 (n + 3) α P N (n). In a potential experimental protocol in which the CM diffusion is monitored by starting from a polymer of a given size n 0 + 3, nonlinear corrections in time mark the crossover to the final Brownian regime in Eq. (17). At variance, in an arrangement in which the initial polymer sizes are distributed according to equilibrium, P N (n 0 , 0) = P N (n 0 ), Eq. (17) turns out to be valid at all time t ≥ 0. Fig. 1 summarizes these behaviors.
It is common practice to quantify the non-Gaussian behavior in terms of the kurtosis (κ X = 3 for Gaussian variables). Note that for simplicity here we refer only to the x-component of the random vector R CM ≡ (X CM , Y CM , Z CM ). From the subordination formula, Eq (15), we get and using again Eq. (12), where we have used the Markov property P N (n , t |n , t , n 0 ) = P N (n , t |n , t ) for t ≥ t ≥ 0. Integrals in Eqs. (16), (20) can be performed, e.g., through Eq. (B.7) in the Appendix, so that the function κ X (t|n 0 ) can be calculated exactly.
Summarizing, we can depict two different scenarios for the CM dynamics of a polymer undergoing polymerization/depolymerization when in contact with a chemostatted monomer bath: For a protocol with an initial specified polymer size n 0 + n min (P N (n 0 , 0) = δ n 0 ,n 0 ), Eqs. (19), (20), (C.1) imply an initial Gaussian behavior with κ X (0 + |n 0 ) = 3. As the distribution in polymer sizes spreads, κ X grows, reaches a maximum, and then returns to κ X (t|n 0 ) ∼ 3 for t τ as specified below. If one starts instead with the equilibrium initial conditions P N (n 0 , 0) = P N (n 0 ), Eq. (20) simplifies which together with Eqs. (19), (C.1) provide the initial kurtosis where Φ(g, α, 3) ≡ ∞ n=0 g n /(n + 3) α is the Lerch transcendent function [46]. Fig. 2 plots Eq. (22) showing that at the critical point g = 1 the equilibrium initial kurtosis diverges for all models of polymer dynamics. Note that for polymer reptation where α = 2, we have a power law divergence (κ X (0 + ) ∼ (1 − g) −1 ), which is logarithmic corrected (κ X (0 + ) ∼ (1 − g) −1 | ln(1 − g)| −2 for the Rouse dynamics (α = 1). In the Zimm model (α = 1/2) the divergence is logarithmic (κ X (0 + ) ∼ − ln(1 − g)) . At large time t τ , the kurtosis does not depend on the initial conditions anymore; through Eq. (21) and (B.7) it is easy to prove a universal power-law decay for the excess kurtosis, independent of the polymer model: κ X (t|n 0 ) − 3 ∝ t τ 1/t. In Fig. 3 we report the evolution of the kurtosis, for both equilibrium and nonequilibrium initial polymer sizes. We observe that close to the critical point, where P N (n) tends to be uniform over an infinite support, the equilibrium results become largely independent of the specific choice of n min . On the contrary, the value of n min affects the nonequilibrium behavior e.g. in situations where a specific initial size n 0 is considered.
Under equilibrium conditions, it is interesting to look at the shape of the initial non-Gaussian PDF for the polymer CM. In order to do so, it is convenient to switch to the unit-variance dimensionless variable X CM (t) ≡ X CM (t)/ E[X 2 (t)]. From Eq. (15), as t → 0 + , we have which only depends on g and α. As displayed in Fig. 4, the tails of this PDF increase with g. At large |x| the PDF is asymptotic to the Gaussian cutoff ∼ e −3 α Davx 2 /(2D 0 ) , and as g → 1 this cutoff is pushed towards |x| → ∞. Consistently with the divergent behavior of the initial kurtosis, tails of the PDF similarly increase with α (see Fig. D3 in the Appendix).
For Zimm and Reptation dynamics, one should also consider intrinsic fluctuations of the diffusion coefficient due to long-range effects triggered by hydrodynamics and entanglement, respectively. For equilibrium initial conditions such effects add a correction [23] to the initial kurtosis which becomes negligible at the critical point g → 1, where the fluctuations of the polymer size due to the polymerization/degradation processes dominate the dynamics. On the contrary, for out-of-equilibrium initial condition the intrinsic fluctuations addressed in Ref. [23] produce a small-time correction to the plots in Fig. 3. Eq. i.c. n 0 =0 n 0 =1 n 0 =10 n 0 =20 asymptotic g=0.96, α=1

Conclusions
We have proposed an analytically solvable microscopic model underpinning Brownian yet non-Gaussian diffusion, a phenomenon ubiquitous in many soft matter and biological systems. The model describes the diffusion dynamics of the CM of a polymerizing/depolymerizing chain in contact with a chemostatted monomer bath. A situation of potential biological importance for such a scenario is that of polymers that bind at specific receptors on cell membranes. In this case the size fluctuations of the tagged particle not only strongly affect the diffusion coefficient of the polymer and thus its mean first passage time to the target, but, once bounded, it may also influence the efficiency of the polymer to trespass the membrane. We have shown that monomer concentration emerges as a natural experimental parameter to study the system under critical conditions. Likewise critical opalescence, we have found that the initial kurtosis κ of the CM location becomes a power-law diverging dynamical response when the monomer concentration c → k − /k + (g → 1). Another conceivable experimental parameter is the concentration of chemically inert macromolecules, added in solution.
Starting at low concentrations, by increasing the crowdness of the environment, it is possible to convert dynamical conditions from the Zimm-hydrodynamic regime to Rouse-diffusion and up to reptation [32,36], thus effectively modifying α. This is yet another possible experimental route to emphasize non-Gaussian behavior. Finally, the established connection with queuing theory may reveal to be useful also where more complex polymerization/depolymerization processes lead to intricate polymer topologies like branched polymers. Extensions of the theory along these directions would be very interesting.

Appendix B. Polymerization -Queuing M/M/1 model
We briefly review the mean-field solution [39] for the P N (n, t|n 0 ) fulfilling ∂ t P N (n, t|n 0 ) = µ P N (n + 1, t|n 0 ) + λ P N (n − 1, t|n 0 ) −(µ + λ)P N (n, t|n 0 ) (n > 0) The probability generating function G(z, t|n 0 ) ≡ ∞ n=0 z n P N (n, t|n 0 ) satisfies U(n) being the (discrete) unit step function and g ≡ λ/µ. Applying the Laplace transform, (·) ≡ ∞ 0 dt e −θ t (·), to Eq. (B.2), one gets Analyzing the zeroes z 1 (θ), z 2 (θ) of the denominator of the last expression and using Rouché's theorem leads to where z 1 is the zero with |z 1 | < 1. At this point G N (z, θ|n 0 ) can be explicitly written as a series expansion, in order to identify P N (n, θ|n 0 ). Finally, application of the inverse Laplace transform provides the desired solution: where I n (z) is the modified Bessel function of the first kind. Using the asymptotic behavior of the modified Bessel function, it can be seen that lim t→∞ P N (n, t|n 0 ) = P N (n). An integral representation of this solution is given by [47] P N (n, t|n 0 ) = U(n) δ n,n 0 − U(n) U(n 0 ) g n−n 0 2 where we point out the identity, obtained by taking the limit t → ∞: √ g cos(θ) .
(B.8) From this expression, the auto-correlation coefficient turns out to be The latter approximation is obtained imposing we obtain the following leading terms for short time evolution of the squared displacement of the polymer CM .
Once an ensemble of simulations has been generated, probability density functions (PDFs) and moments of the stochastic process R CM (t) can then be numerically inferred. The following plots summarize comparisons between numerical and analytical results with dt/τ = 0.1.  Eq. i.c-Analytic Eq. i.c-Simulation n 0 =10-Analytic n 0 =10-Simulation g=0.96, α=1 Figure D2. Time evolution of the CM x-kurtosis obtained through Gillespie-Langevin simulations (symbols) contrasted with analytical results for both equilibrium and nonequilibrium initial polymer sizes. Simulation data are obtained averaging over 9 × 10 5 (n 0 = 10) and 4 × 10 5 realizations (equilibrium).  Figure D3. Unit-variance initial x-PDF for the CM of different polymer models with equilibrium sizes P N (n 0 , 0) = P N (n 0 ). Gillespie-Langevin simulations (symbols) are contrasted with analytical results. Simulation data are obtained averaging over 4 × 10 5 realizations. For comparison purposes, a unit-variance Gaussian PDF is also plotted in red dotted line.