Time and space generalized diffusion equation on graph/networks

Time and space generalized diffusion equation on graph/networks Fernando Diaz-Diaz and Ernesto Estrada arXiv:2202.00318v1 [physics.soc-ph] 1 Feb 2022 Institute of Cross-Disciplinary Physics and Complex Systems, IFISC (UIB-CSIC), 07122 Palma de Mallorca, Spain Normal and anomalous diffusion are ubiquitous in many complex systems [1]. Here, we define a time and space generalized diffusion equation (GDE), which uses fractional-time derivatives and transformed d-path Laplacian operators on graphs/networks. We find analytically the solution of this equation and prove that it covers the regimes of normal, sub- and superdiffusion as a function of the two parameters of the model. We extend the GDE to consider a system with temporal alternancy of normal and anomalous diffusion which can be observed for instance in the diffusion of proteins along a DNA chain. We perform computational experiments on a one-dimensional system emulating a linear DNA chain. It is shown that a subdiffusive-superdiffusive alternant regime allows the diffusive particle to explore more slowly small regions of the chain with a faster global exploration, than a subdiffusive-subdiffusive regime. Therefore, an alternancy of sliding (subdiffusive) with hopping and intersegmental transfer (superdiffusive) mechanisms show important advances for protein-DNA interactions. 1. INTRODUCTION Diffusion–the net movement of particles in an environment, generally from a region of higher concentration to a region of lower concentration–is ubiquitous in natural and man-made systems. In the absence of obstacles to diffusion and traps, the diffusive particles describe a random walk motion on the environment, such that their mean squared displacements (MSD) scale linearly with time, hx2 i ∝ t. This process is known as normal diffusion. However, the existence of obstacles in the environment may trigger long-jumps of the diffusive particle [2–4], such that the mean displacement of the particles is bigger than that of the normally diffusing ones in the same period of time, i.e., hx2 i ∝ tγ>1 . This type of process is known as superdiffusive. On the other hand, it is possible that the environment has regions acting as traps for the particles, where they are retained for longer times than in a normal diffusive process. In this case, hx2 i ∝ tγ<1 , and the process is known as subdiffusive. From the modeling perspective there are several approaches to describe anomalous (sub- and super-) diffusion [5, 6]. From a physical perspective these processes have analogues in terms of anomalous heat conduction [7], where normal diffusion implies normal heat conduction, superdiffusion implies anomalous heat conduction with a divergent thermal conductivity and subdiffusion implies anomalous heat conduction with a convergent thermal conductivity. In the last case the system is a thermal insulator in the thermodynamic limit. It has been remarked that subdiffusion may arise as the result of the coexistence of time-periods dominated by normal transport with periods in which there is no effective transport. The last can emerge when the diffusive particle is temporarily trapped as a result of geometrical complexity and interactions with the environment. This could be clearly the case of the travel of contaminants in groundwater, which display much longer times than the ones expected from the classic diffusion. The motion of proteins while sliding on DNA during target search is believed to be subdiffusive in general [8]. Simulations results for the case of T7 RNA polymerase promoter search on T7 DNA has been found to be subdiffusive for short times and asymptotically approaching normal diffusion [9]. More recent intensive computational simulations also pointed out the important role of subdiffusive process in the diffusive search of proteins for their specific binding sites on DNA in the presence of the macromolecular crowding in cells [10]. Nowadays it is well established that the molecular crowding of the internal cellular environment induces the emergence of anomalous subdiffusion of cytoplasmic macromolecules. This has been verified by means of fluorescence correlation spectroscopy and computer simulations [11], fluorescence correlation spectroscopy [12], and by tracking fluorescently labeled mRNA molecules [13]. The complexities of the process where studied with globular proteins dispersed in aqueous solution of poly(ethylene oxide) (PEO) to mimic a crowded environment. Using state-of-the-art neutron spin echo (NSE) and small-angle neutron scattering (SANS) techniques it was observed a fast dynamic corresponding to diffusion inside a trap built by the polymer mesh with slower process corresponding to the long time diffusion on macroscopic length scales [14]. It has also been found that water molecules jump randomly between trapping sites on protein surfaces, giving rise to subdiffusion. At longer times the subdiffusive exponent gradually increases towards normal diffusion due to a many-body volume-exclusion effect [15]. The intermittency of fast and slow processes as the one described by the globular proteins in a crowded environment [14] can also be found in other scenarios. For instance, it has been reported that a proliferating, diffusing tumor within different surrounding tissue conditions migrates not only by using normal diffusion, but also using combinations of subdiffusion, superdiffusion, and even ballistic diffusion, with increasing mobility of the tumor cell when haptotaxis 2 and chemotaxis toward the host tissue surrounding the proliferative tumor are involved [16]. The cytoskeleton (CSK), a crowded network of structural proteins that stabilizes cell shape and drives cell motions, displays spontaneous subdiffusive bead motions at short times followed by superdiffusive motion at longer times. The intermittency of the motions depended on both the approach to kinetic arrest and energy release due to ATP hydrolysis [17]. If we consider the protein’s diffusive transport on DNA in a wider perspective, i.e., not only considering the sliding process, then we observe a whole range of normal and anomalous behavior. Using recent advances in single molecule detection it has been observed that proteins diffusing along DNA follow different mechanisms, such as (i) random collision, (ii) sliding, (iii) hopping, (iv) intersegmental transfer and (v) active translocation [18, 19]. While the sliding can give rise to subdiffusive and normal diffusion behavior, the hopping is known in other systems to produce superdiffusive behavior [20–23]. Additionally, intersegmental transfers can transport a protein from one site in the DNA to another very distant from the original one [24–26], which can give rise to superdiffusive behavior. Many of the complex systems in which these normal and anomalous diffusion processes take place form interaction networks [27, 28]. For instance, DNA can be represented as a linear chain on which a protein is diffusing between its nodes and edges. Therefore, here we consider the time and space generalization of the diffusion equation on graphs/networks. We consider time-fractional derivatives, which account for nonlocality by time or dynamic memory [29, 30], and long-range jumps in the graph/network through the d-path Laplacian operators [31–33]. We first define the generalized diffusion equation (GDE), study the main properties of its solution and find analytically the conditions for the existence of normal, sub- and superdiffusion. We also consider a time-varying GDE such that the three diffusive regimes appear intermittently with time. Finally, we apply this approach to the study of the diffusion of a particle through a linear chain representing a protein diffusing through DNA. 2. PRELIMINARIES Here we will use interchangeably the terms graph and network for G = (V, E), which in general will be connected and undirected. If the number of nodes n = #V is infinite we then assume that G is locally finite (all vertices have finite degree). Furthermore, let l2 (V ) be the Hilbert space of square-summable functions on V . We begin by defining the d-path Laplacian operator on G. Let dist(v, w) be the length of the shortest path between v and w, and let kd (v) be the d-path degree of the vertex v, defined by: kd (v) = #{w ∈ V : dist(v, w) = d}. (2.1) Let f be a function acting over the set of vertices of G. Then, the d-path Laplacian operator Ld is defined by: X (Ld f )(v) = (f (v) − f (w)). (2.2) w∈V :dist(v,w)=d Let ev (w) be the orthonormal basis ev (w) = ( 1 if w = v, 0 otherwise. (2.3) Then, we have   kd (v) (Ld ev )(w) = −1  0 if w = v, if dist(v, w) = d, otherwise. (2.4) The Mellin-transformed d-path Laplacian is the following weighed sum of d-path Laplacians: L̃(s) = ∆ X Ld d−s , (2.5) d=1 where ∆ is the diameter of the graph. This operator preserves several key properties of the graph Laplacian [31, 32]: Lemma 1. The Mellin-transformed d-path Laplacian is positive semidefinite. Furthermore, let kd,max := max{kd (v) : v ∈ V }. If kd,max ≤ Ck α , then L̃(s) is bounded for all s ∈ C with R(s) > α + 1. 3 When s → ∞, all non-unity entries of L̃(s) vanish and we have the standard graph Laplacian operator L. The physical difference of the d-path Laplacian operator with the fractional powers of the standard Laplacian have been analyzed in [34]. The standard diffusion equation on a graph is formulated on the basis of this operator as follows: ∂f (t) = −DLf (t),f (0) = f0 , ∂t (2.6) where D is the diffusion coefficient, hereafter taken always to be unity. This equation has the well-known solution: f (t) = e−Lt f0 . (2.7) The function f (t) can be interpreted as a probability density function (pdf) of the position of a hypothetical diffusing particle. Indeed, the diffusion equation ensures that several key properties of a pdf are satisfied: (i) If (f0 )n ≥ 0 ∀n, the exponential operator P e−Lt ensures that the time evolved pdf remains positive: (f (t))n ≥ 0 ∀n, t. (ii) If the initial condition is normalized, ( n (f0 )nP= 1) then f (t) remains normalized for any t. This is because the diffusion equation preserves the 1-norm k f (t) k1 := n fn (t): ! X X X X X d k f (t) k1 d X Ljj fj − Lij fj = Lij fj = 0, (2.8) (Lf )i = − Lij fj − fi = − = dt dt j i i i,j6=i i,j6=i i,j6=i P where we have used Ljj = − i6=j Lij . The above equation implies that k fi k1 (t) is constant, as previously stated. When the diffusion equation is defined as before, the mean square displacement (MSD) hx2 i of the diffusive particle scales linearly with time: hx2 i ∝ t, independently of the initial conditions and the diffusion coefficient. This is known as normal diffusion in the literature. However, many physical systems display the so-called anomalous diffusion [5, 6], where the MSD scales as a power law with time: hx2 i ∝ tγ . If γ < 1, the system is in a subdiffusive regime, whereas if γ > 1, the system is in a superdiffusive regime. Moreover, in a normal diffusive regime the maximum of the probability density function (pdf), fmax , decreases as fmax (t) ∝ t−0.5 , while in superdiffusive and subdiffusive ones it decays as fmax (t) ∝ t−γ>0.5 and fmax (t) ∝ t−γ<0.5 , respectively. Similar decays exist also for the so-called Full Width at Half Maximum (FWHM) of the pdf, namely FWHM ∝ tγ>0.5 and FWHM ∝ tγ<0.5 for super- and subdiffusive regimes. 3. TIME AND SPACE GENERALIZED DIFFUSION ON GRAPHS/NETWORKS Here we define a time and space generalized diffusion equation on graphs/networks in the following way. Let Dtα be the Caputo fractional derivative and let L̃(s) be the Mellin-transformed d-path Laplacian operator on G. Then, the generalized diffusion equation (GDE) is defined as Dtα f (t) = −L̃(s)f (t),f (0) = f0 , (3.1) where Dtα f (t) 1 = Γ(1 − α) Z t 0 f ′ (τ ) dτ , (t − τ )α (3.2) f ′ (τ ) denotes the usual derivative. Here, 0 < α ≤ 1 and 0 < s < ∞. Obviously, when s → ∞, we have Dtα f (t) = −Lf (t), where L is the standard graph Laplacian. This equation accounts for a time-fractional process only, without any spatial long-range jumps. On the other hand, when α = 1 we get dfdt(t) = −L̃(s)f (t), which accounts for long-range spatial jumps in the graph as studied previously in [32]. Our first result is the general expression for the solution of the GDE defined before. We state this result in the following. Theorem 2. The solution of the GDE (3.1) is given by f (t) = Eα (−L̃(s)tα )f0 , (3.3) where Eα (...) is the Mittag-Leffler function of the corresponding matrix: Eα (−L̃(s)tα ) = ∞ X (−L̃(s)tα )j j=0 Γ(αj + 1) . (3.4) 4 Proof. Let L̃(s) = U ΛU † , where Λ = diag(λ1 , ..., λn ) and U is a unitary matrix of eigenvectors of L̃(s). The dagger symbol denotes the conjugate transpose. Thus: Dtα f (t) = −U ΛU † f (t). (3.5) Defining y(t) = U † f (t), we obtain a set of decoupled fractional differential equations: Dtα yi (t) = −λi yi (t), † yi (0) = (U f0 )i =: y0i . (3.6) (3.7) Let us use the Laplace transform of this fractional differential equation, and using the properties L{Dtα y}(u) = α−1 uα L{y}(u) − uα−1 y(0) as well as L{Eα (−atα )}(u) = uuα +a , we find the solution: yi (t) = y0i Eα (−λi tα ). (3.8) Finally, we undo the change of basis y(t) = U † f (t) to obtain: f (t) = Eα (−L̃(s)tα )f0 , (3.9) where the definition of a matrix function g acting on a diagonalizable matrix A = U ΛU † is: g(A) = U g(Λ)U † , and g(Λ) := diag(g(λi )). A. Analysis of different diffusive regimes We now focus on understanding how the GDE can give rise to different diffusive regimes in a given graph. In particular, we will focus here only in the one-dimensional case. For that, we will consider an infinite path graph, P∞ , which corresponds to the one-dimensional case. Anomalous diffusion in one-dimension has received much attention in the literature [35–39] due to its practical relevance. Let us start by defining the Fourier transform operator F and its inverse F −1 , respectively: 1 X F{f }(k) = √ fn eink , (3.10) 2π n∈Z Z π 1 e−imk g(k)dk. (3.11) (F −1 {g})m = √ 2π −π Let us first obtain the solution of the GDE for the P∞ graph. Lemma 3. The fundamental solution of the GDE for the infinite path graph P∞ and initial condition fn (0) = f0 = δn0 , can be expressed as: Z π 1 fm (t) = dke−ikm Eα (−tα ls (k)), (3.12) 2π −π with where ζ(s) = P∞ 1 k=1 ks ls (k) = 2ζ(s) − Lis (eik ) − Lis (e−ik ), P∞ n is the Riemann zeta function and Lis (x) = n=1 xns is the polylogarithm. (3.13) Proof. The Fourier transform of the Mellin-transformed d-path Laplacian operator L̃(s) is [32]: F{L̃(s)fn } = ls (k)F{fn }, (3.14) F{f (t)} = F{Eα (−L̃(s)tα )f0 } = Eα (−ls (k)tα )F{f0 }. (3.15) therefore: Using the inverse Fourier transform and substituting the initial condition fn (0) = (f0 )n = δn0 , we finally arrive at the expression: Z π 1 fm (t) = dke−ikm Eα (−tα ls (k)). (3.16) 2π −π 5 The time scaling of this expression is unclear, due to the polylogarithms inside the Mittag-Leffler function. However, to determine the diffusive regime that the system exhibits, it is sufficient to characterize the behavior of the solution in the limit of long times. Thus, we perform an asymptotic approximation, which is accounted for in the following result. Lemma 4. Let β > 0 and let l : [−π, π] → R be a continuous function satisfying: l(k) > 0 f or k ∈ [−π, π]\0 l(k) ∼ c|k| β (3.17) as k → 0 (3.18) A(k) with some c > 0. A(k) ∼ B(k) denotes that limk→0 B(k) = 1. Then:   Z π Z π 1 1 1 −ikm α −ikm α β −1 α β √ Eα (−ct |k| ) e Eα (−t l(k))dk → e Eα (−ct |k| )dk = F 2π −π 2π −π 2π (3.19) as t → ∞. Proof. Let ǫ > 0 and let us decompose the Fourier transform integral into: Z π Z ǫ Z 1 1 1 e−ikm Eα (−tα l(k))dk = e−ikm Eα (−tα l(k))dk + e−ikm Eα (−tα l(k))dk. 2π −π 2π −ǫ 2π [−π,−ǫ]∪[ǫ,π] Let us take the limit ǫ → 0, so that the first integral becomes: Z ǫ Z ǫ e−ikm Eα (−tα l(k))dk → e−ikm Eα (−tα c|k|β )dk. −ǫ (3.20) (3.21) −ǫ The second integral from (3.20) is negligible: Z [−π,−ǫ]∪[ǫ,π] e−ikm Eα (−tα l(k))dk ≤ Z Eα (−tα l(k))dk = 2 [−π,−ǫ]∪[ǫ,π] Z π ǫ Eα (−tα l(k))dk → 0 as t → ∞. (3.22) In particular, for a given ǫ, the integral is negligible whenever t ≫ l(ǫ)−1/α . For identical reasons, the integral Z e−ikm Eα (−tα c|k|β )dk (3.23) [−π,−ǫ]∪[ǫ,π] is also negligible when t → ∞. Thus, we can express the Fourier transform integral as:   Z π Z π 1 1 1 e−ikm Eα (−tα l(k))dk → e−ikm Eα (−tα c|k|β )dk = F −1 √ Eα (−tα c|k|β ) . 2π −π 2π −π 2π Corollary 5. The solution of the GDE, when t → ∞, is given by:   Z π 1 1 α β −ikm α β −1 √ Eα (−ct |k| ) , e Eα (−ct |k| )dk = F fm (t) = 2π −π 2π (3.24) (3.25) with β= ( s−1 2 if 1 < s < 3, if s > 3. (3.26) Proof. According to Theorem 6.5 of [32], the function ls (k) = 2ζ(s)−Lis (eik )−Lis (e−ik ) has the following asymptotics as k → 0: ( π − Γ(s) cos |k|s−1 if 1 < s < 3, ( sπ 2 ) ls (k) ∼ (3.27) ζ(s − 2)k 2 if s > 3, Moreover, ls (k) is continuous and fulfills ls (k) > 0 ∀k 6= 0, ls (0) = 0 ([32], Lemma 5.1). Thus, ls (k) fulfills the hypotheses of Lemma 4, so its asymptotic approximation is given by eq. (3.19). 6 Equation (3.25) fully determines the solution of the GDE for long times. Unfortunately, the Fourier transform cannot be solved in terms of elementary functions. Nevertheless, the time scaling can be obtained with the following Theorem. α Theorem 6. The solution of the GDE f (t) has the asymptotic scaling with time f (t) ∝ t− β , as t → ∞. Proof. Let us exploit the parity of the integrand in eq. (3.25) and expand the cosine into a power series, we obtain: Z X1Z π 1 π (km)2j fm (t) = Eα (−ck β tα ). (3.28) dk cos(km)Eα (−ck β tα ) = dk(−1)j π 0 π (2j)! 0 j We perform the change of variables u = ck β tα and split the integration interval into two: Z ∞  Z ∞ Z π 2j+1−β 2j+1−β 1 2j β α β β duu duu Eα (−u) − Eα (−u) . dkk Eα (−ck t ) = 2j+1 0 π β ctα 0 β(ctα ) β (3.29) π β ctα The second integral is bounded by the constant e− 2 , as we will show now. For sufficiently big u, Eα (−u) < e−u 2j+1−β u and u β < e 2 , so: Z ∞ Z ∞ 2j+1−β u π β ctα Eα (−u)u β du < (3.30) e−u e 2 du = e− 2 . π β ctα π β ctα β α − π 2ct approaches zero and thus the second integral from eq. (3.29) is negligible; therefore:   X (−1)j m2j 1 I0 (α, β) − αβ m2 I1 (α, β) − 3α 1 β + ... fm (t) = t t I (α, β) = − , (3.31) j 1 3 π (2j)! β(ctα ) 2j+1 πβ 2 β cβ cβ j In the limit t → ∞, e R∞ 2j+1−β where Ij (α, β) = 0 duu β Eα (−u). For large times and small m (i.e., in the central region), the dominant term in the sum is the one with the smallest absolute value of the exponent of t; i.e., the term j = 0. Hence, the asymptotic scaling with time is α fm (t) ∝ t− β . (3.32) Additionally, let us find the scaling of two relevant observables: the height of the maximum of the pdf and the Full Width at Half Maximum (FWHM). Lemma 7. Theα maximum of the pdf, fmax , and FWHM fulfill the following asymptotic scaling laws in time: i) fmax (t) ∝ t− β , α ii) FWHM(t) ∝ t β Proof. The proof of (i) is straightforward using eq. (3.31) for m = 0: fmax (t) = α 1 I0 (α, β) − αβ ∝ t− β . t 1 πβ c β (3.33) α To prove (ii), we make the rescaling m → mt β in eq. (3.25). The function f (m, t) becomes: Z π Z π α α α 1 1 −α −ikt β m α β β β e e−ik̃m Eα (−c|k̃|β )dk̃ = t− β f (m, 1), Eα (−ct |k| )dk = t f (mt , t) = 2π −π 2π −π (3.34) α where k̃ = t β k. Let us apply this equation to find the scaling of the FWHM. Let ξ(t) be the position at which the pdf reaches its Half Maximum at time t (i.e., the ”Half Width at Half Maximum”), so that FWHM(t) = 2ξ(t). The (t) FWHM can be determined through the equation f (ξ(t), t) = fmax . Using the time scaling of the maximum obtained 2 in (i), we can express this as: α α fmax (1)t− β f (ξ(t), t) = = f (ξ(1), 1)t− β . 2 (3.35) α Comparing this equation with the scaling relation (3.34), it follows that ξ(t) = ξ(1)t β , and thus: α FWHM ∝ t β . (3.36) 7 FIG. 3.1: Phase diagram with the three possible diffusive regimes (superdiffusion, subdiffusion and standard diffusion) as a function of the parameters α and s. Remark: The MSD for processes with infinite variance of the step size distribution diverges [40, 41], unless the step size distribution and the waiting time distribution are correlated [42]. This is a well-known limitation of Levy flights and other stochastic processes that lead to superdiffusion [41]. Proposals of pseudo mean square displacements that remain finite for superdiffusive processes have been made [5, 40], but their analytical calculation is far from trivial. 1 1 Let us now recall that for normal diffusion, the corresponding scaling laws are: i) f (t) ∝ t− 2 , ii) fmax (t) ∝ t− 2 ,and 1 iii) FWHM(t) ∝ t 2 . Let us focus here on the case 1 < s < 3, where β = s − 1. We find that the three scaling laws predict the same location of the three possible diffusive regimes: Standard diffusion: α = Superdiffusion: α > Subdiffusion: α < s−1 2 . s−1 2 . s−1 2 . For s > 3, β = 2, which implies that subdiffusion happens for any α ∈ (0, 1), and standard diffusion for α = 1. A phase diagram with the different diffusive regimes can be found in Figure 3.1. 4. TIME-DEPENDENT GENERALIZED DIFFUSIVE PROCESSES Let us consider here a system in which the different diffusive regimes alternate with time. That is, a system in which the diffusive particle behaves subdiffusively at a given time, then alternates with normal diffusion and superdiffusion as the time goes on. These processes can alternate cyclically as we will see later on this work. 8 A. Stepwise temporal dependence The simplest way of achieving this is to let the parameter s depend on time. This idea was previously explored by Allen-Perkins et al. [43] for the case of the d-path Laplacian operator. Let us define: X L̃(s(t)) = Ld d−s(t) , (4.1) ( s1 if t < tsw , (4.2) s(t) = s2 if t > tsw . We have the following generalized diffusion equation: Dtα f (t) = −L̃(s(t))f (t) , f (0) = f0 . (4.3) For a step-like s like eq. (4.2), one can split the system into two time-independent equations: Dtα fA (t) = −L̃(s1 )fA (t) Dtα fB (t) = −L̃(s2 )fB (t) if t < tsw , (4.4) if t > tsw , (4.5) with initial condition fA (0) = f0 . Then, expressing L̃(sq ) (q ∈ {A, B}) as L̃(sq ) = Uq Λq Uq† , where Uq are unitary and Λq = diag((λq )i ), and defining yq = Uq† fq , we get a set of decoupled differential equations. The solution of each one is: (yq )i (t) = Eα (−(λq )i tα )(C̃q )i , (4.6) where C̃q are integration constants. C̃A can be obtained using the initial condition fA (0) = f0 and C̃B is found imposing the continuity of the function at t = tsw : fA (tsw ) = fB (tsw ). Undoing the change of basis induced by Uq† , we obtain the following solution in terms of matrix functions: ( Eα (−L̃(s1 )tα )f0 if t ≤ tsw , f (t) = (4.7) α α −1 α Eα (−L̃(s2 )t )(Eα (−L̃(s2 )tsw )) Eα (−L̃(s1 )tsw )f0 if t > tsw . From now on, we assume that t > tsw . We express eq. (4.7) as: f (t) = Eα (−L̃(s2 )tα )f1 , (4.8) −1 where f1 (tsw ; s1 , s2 ) = Eα (−L̃(s2 )tα Eα (−L̃(s1 )tα sw ) sw )f0 . The introduction of f1 reduces the problem to a timeindependent expression, analogous to eq. (3.9), but with a different initial condition. The time-evolved system is then Z  1 X π fm (t) = dke−ik(m−n) Eα (−tα ls2 (k))(f1 )n → F −1 Eα (−c2 tα |k|β2 )F{f1 }(k) (4.9) 2π n −π as t → ∞, where we have used the asymptotic equivalence ls2 (k) ∼ c2 |k|β2 as k → 0 and Lemma 4. The time scaling of the pdf can be calculated using the procedure from Theorem 6, obtaining: fm (t) = X (−1)j (m − n)2j 1 2j+1 Ij (α, β2 )(f1 )n , π (2j)! β2 (c2 tα ) β2 j,n (4.10) α which confirms that the dominant time scaling is still t− β2 . In other words, the phase diagram from the previous section remains unchanged: the system behaves asymptotically as if the Laplacian L(s1 ) had never affected the dynamics. B. Periodic temporal alternancy Consider the following time-dependent parameter s, with period T : ( s1 if nT ≤ t < (n + 21 )T, s(t) = s2 if (n + 12 )T ≤ t < (n + 1)T, (4.11) 9 where n ∈ N ∪ {0}. In order to obtain the full solution of the problem, let us define f (n) (t) as the restriction of f (t) to t ∈ [nT, (n + 1)T ): f (n) (t) = f (t)|[nT,(n+1)T ) . To solve the problem, it suffices to find an expression for f (n) (t) for arbitrary n. Theorem 8. The solution of the GDE with a periodic parameter s(t) and initial condition f (0) = g0 is given by: ( Eα (−L(s1 )tα )gn if nT ≤ t < (n + 21 )T

f (n) (t) = α α T T if (n + 12 )T ≤ t < (n + 1)T , Eα (−L(s2 )tα )(Eα (−L(s2 ) nT + 2 ))−1 Eα (−L(s1 ) nT + 2 )gn (4.12) where: α α   T T ))−1 Eα (−L(s1 ) nT − )gn−1 = gn :=Eα (−L(s1 )(nT )α )−1 Eα (−L(s2 )(nT )α )(Eα (−L(s2 ) nT − 2 2   α α n Y T T −1 α −1 α = )) Eα (−L(s1 ) jT − )g0 . (4.13) Eα (−L(s1 )(jT ) ) Eα (−L(s2 )(jT ) )(Eα (−L(s2 ) jT − 2 2 j=1 Proof. f (0) (t) is given by eq. (4.7), for tsw = T2 . The remaining f (n) (t) can be obtained by induction. Suppose f (n) (t) is given by eq. (4.12). The solution to the GDE in the interval [(n + 1)T, (n + 2)T ) is: ( Eα (−L(s1 )tα )C1 if (n + 1)T ≤ t < ((n + 23 )T , (n+1) f (t) = (4.14) Eα (−L(s2 )tα )C2 if (n + 32 )T ≤ t < (n + 2)T . where C1 , C2 are integration constants. Let us impose the continuity of the function at t = (n + 1)T and t = (n + 32 )T to obtain the values of the integration constants: C1 =Eα (−L(s1 )((n + 1)T )α )−1 Eα (−L(s2 )((n + 1)T )α )×  α α  T T × (Eα (−L(s2 ) nT + ))−1 Eα (−L(s1 ) nT + )gn = gn+1 , 2 2 3 3 C2 =(Eα (−L(s2 )((n + )T )α )−1 Eα (−L(s1 )((n + )T )α )gn+1 . 2 2 (4.15) (4.16) This corresponds to eq. (4.12) for n → n + 1. Thus, by induction, expression (4.12) holds for any n.
α The pdf given by (4.10) is still valid, with f1 = gn if nT ≤ t < ((n + 12 )T and f1 = (Eα (−L(s2 ) nT + T2 ))−1 ×
α ×Eα (−L(s1 ) nT + T2 )gn if (n + 12 )T ≤ t < (n + 1)T . This implies that the time scaling at time t depends only on the parameter s of the system at that time. Consequently, the phase diagram of fig. 3.1 is valid, adapting s to the Mellin transform parameter that acts at the desired time. 5. GENERALIZED DIFFUSION ON A LINEAR CHAIN We consider here the diffusion of a particle along a linear chain or path graph Pn . The motivation of these simulations is provided by the diffusion of a protein along a DNA chain. We consider the double chains of DNA (see Fig. 5.1 a) as the path graph Pn . In doing so, we consider a pair of DNA bases, one from the 5’ chain and another from the 3’ chain, as a node of the path graph. The protein is then modeled as a particle diffusing on the linear chain. We consider here that the protein can diffuse by a combination of the following mechanisms: (i) one-dimensional diffusion or sliding, involving a random walk along the DNA without dissociation as illustrated in Fig. 5.1 b; (ii) jumping, where a protein moves over longer distances via dissociation and then rebinding at a distal location (see Fig. 5.1 c); (iii) intersegmental transfer, involving movement from one site to another via a looped intermediate (Fig. 5.1 d). The first process gives rise to normal diffusion, but we consider that at certain DNA regions the protein makes a longer exploration/repair which acts as a trap producing subdiffusive motion, which can be modeled by varying the exponent α of the Caputo fractional derivative. The mechanisms (ii) and (iii) clearly represent long-jumps controlled 10 FIG. 5.1: (a) Representation of DNA chains and the transcription factor (TF). In the following panels DNA is represented as a chain graph where every pair of bases is represented by a node and two pairs of consecutive nodes are connected by an edge. (b-d) Different processes by which TF inspect and repair DNA: (b) sliding, (c) hopping and (d) intersegmental transfer. by the Mellin exponent in the transformed d-path Laplacian. We should remark that other approaches have been reported in the literature to account for some aspects of the protein diffusion on DNA [44–46]. For the simulations we consider a path graph with N = 1001 nodes. We always take as initial condition fn (0) = δn,501 , i.e., we locate all the diffusive particle at the center of the linear chain. To obtain the probability density functions, we calculated eq. (3.3) using the algorithm from [47] for matrix Mittag-Leffler functions. We then fitted the time evolution of the maximum of the pdf and FWHM to a power law tγ , and compared the numerical exponent γ with the corresponding theoretical prediction. We have noticed that in the superdiffusive regime it is hard to precisely quantify the MSD, since the asymptotic behavior is reached for large times and erratic behaviors can appear in the transient regime (see [48] for a detailed discussion). Moreover, the FWHM method, while providing correct predictions, introduces a higher degree of inaccuracy in the exponent γ. Because of this, we use the method of the decay of the maximum to find the numerical exponent γ. Our first target is to obtain a contour plot indicating how the exponent γ defining the type of diffusive regime changes with the values of the model parameters α and s. The results are illustrated in Fig. 5.2. We recall that for fmax (t) ∝ t−γ : i) γ < 0.5 represents subdiffusion; ii) γ = 0.5 represents normal diffusion; iii) γ > 0.5 represents superdiffusion. Therefore, we can see the curve with the value γ = 0.5 in the Fig. 5.2, which indicates the normal diffusion regime. When we move towards the right-upper corner of the plot, i.e., s → ∞ and α → 1, we are moving to a regime completely dominated by normal diffusion, which is represented by the standard graph diffusion equation. However, it is remarkable that even with strong long-jumps, e.g., s ≈ 2, we can still find regions of normal diffusion if the temporal memory parameter α is relatively small, e.g., α ≈ 0.6. Over this curve we are in a superdiffusive regime and below it we are in the subdiffusive one. It is remarkable that there is an abrupt change of behavior at s = 3, which matches the theoretical prediction that after this point the superdiffusive regime no longer exists. It can be seen that for s < 3 the contour lines of γ abruptly decay, while for s > 3 they are almost parallel to the x-axis. The main conclusion of this result is that, as proved analytically, the time and space GDE accounts for the three different diffusive regimes that can be observed in the diffusion of a protein on DNA: subdiffusion, normal diffusion and superdiffusion. In Fig. 5.3 we give an example of the time evolution of the diffusion of a particle along a linear 11 FIG. 5.2: Exponent γ as a function of the control parameters s and α. γ > 0.5 corresponds to superdiffusion and γ < 0.5, to subdiffusion. The exponents γ were obtained by fitting the decay of the maximum to a power law. For β each pair (s, α), the time range used to measure this decay was t ∈ [0, 10 α ], where β is given by eq. (3.26). chain of N = 1001 nodes. We use α = 0.5, s = 2.5, which according to the contour plot in Fig. 5.2, corresponds to the subdiffusive regime. The value of γ calculated from the slope of the curve in Fig. 5.3 (b) is γ ≈ 0.33 ± 0.01 and according to the one of Fig. 5.3 (c) is γ ≈ 0.31 ± 0.02. The analytical value is γ = 1/3. We now turn our attention to the study of the time-varying GDE with periodic variation of the parameter s. We consider that s oscillates between the values s1 = 200, and s2 = 2. That is, for a given values of α the particle starts diffusing without long-jumps (s1 = 200). Then, after half a period of time it swaps to a regime where long-jumps are allowed (s2 = 2), and the process is repeated cyclically. We first fix α = 0.5 and obtain the results illustrated in Fig. 5.4 (a). For 0 < t < 2 we have fmax (t) ≈ 0.3887t−0.1987 and for 2 < t < 4 we obtain fmax (t) ≈ 0.4381t−0.3982 . That is, when α = 0.5 the diffusive particle oscillates between two subdiffusive regimes. In fact, the global fit of the process is given by: fmax (t) ≈ 0.3671t−0.2472 , which is a clear signature of subdiffusion. We now consider the case where s1 = 200 and s2 = 2 as before, but using α = 0.9 as illustrated in Fig. 5.4 (b). Here, fmax (t) ≈ 0.3368t−0.3513 for 0 < t < 2 and fmax (t) ≈ 0.5108t−1.133 for 2 < t < 4. The global process has scaling fmax (t) ≈ 0.2715t−0.4651 . That is, the global process is a subdiffusive regime, although the diffusive particle alternates between a subdiffusive regime, i.e., γ ≈ −0.3513 for 0 < t < 2 , and a superdiffusive motion, γ ≈ −1.133 for 2 < t < 4. The importance of this difference between the local temporal scale and the global one is revealed when we plot the decay of fmax vs. t for both values of the fractional parameter α as illustrated in Fig. 5.5. As can be seen, the process where α = 0.5 (blue circles) goes initially much faster than that where α = 0.9 (red triangles). Then, there is a time in which the process with α = 0.9 is much faster than that with α = 0.5. In other words, the combination of subdiffusion with superdiffusion, like when α = 0.9, allows the particle to make a slower initial exploration of a region of the linear chain in comparison with a subdiffusive-subdiffusive exploration. Additionally, the subdiffusive-superdiffusive process produces a faster global convergence of the process due to the long-jumps occurring in the superdiffusive regime. Notice that the subdifusive-subdifusive regime is obtained when the fractional parameter is relatively small, which correspond to systems with relatively large temporal memory. However, the subdiffusive-superdiffusive alternancy is obtained when the temporal memory is relatively small. Translating these results to the case of a protein diffusing along a DNA chain they mean that the alternant combination of sliding with jumping and/or intersegmental transfer offer some important advances to the exploration of the DNA by the protein. In this case, the slow subdiffusive regime 12 Probability 0.3 0.2 0.1 0 480 500 520 Node (a) (b) (c) FIG. 5.3: Illustration of the time evolution of a diffusive particle on a path graph with N = 1001 nodes, using α = 0.5, s = 2.5. The pdf is initially a delta distribution at the middle of the linear chain. (a) Probability density function in the different nodes, for three different times: t = 1 (solid line),t = 10 (dashed line) and t = 100 (dotted line). (b) Decay of the height of the maximum with time. (c) Increase of the FWHM as a function of time. allows a detailed exploration of small DNA regions to find potential targets and the fast superdiffusive regimes allow an exploration of vast regions of the DNA chain in relatively short times. 6. CONCLUSIONS AND OUTLOOK We have defined a time and space GDE on undirected graphs/networks. It uses a combination of fractional derivatives and Mellin-transformed d-path Laplacian operator. We have found analytically the solution of this equation and obtained the regions of the parametric space for which an infinite one-dimensional system displays, normal, suband superdiffusion. We have illustrated how this GDE can be applied to the study of the diffusion of proteins along the one-dimensional structure of DNA, where the mechanisms of sliding, hopping and intersegmental transfer, may give rise to normal, sub- and superdiffusive behaviors. We have also considered a GDE in which the parameters of the model change with time allowing the temporal alternancy of the normal and anomalous diffusive regimes. The current model is useful for any discrete system in which any combination of the normal and anomalous diffusive regimes exists. The extension of this model to consider directed graphs, multigraphs and simplicial complexes can be performed to extend the areas of application of this approach. Also important should be the analysis of networks beyond the one-dimensional case presented here and to consider the influence of network topologies, e.g., smallworldness, scale-freeness, etc., on the diffusive dynamics. All in all we consider that this GDE will open new research avenues in the study of dynamical processes on graphs/networks. ACKNOWLEDGEMENT FD-D and EE acknowledge support from the Spanish Agency of Research (AEI) through Maria de Maeztu Program for units of Excellence in R&D (MDM-2017-0711). FD-D thanks finantial support MDM-2017-0711-20-2 funded by MCIN/AEI/10.13039/50110 0 011033 and by FSE invierte en tu futuro. EE thanks Grant PID2019-107603GB-I00 by MCIN/AEI /10.13039/50110 0 011033. [1] A. Bunde, J. Caro, J. Kaerger, and G. Vogl, Diffusive Spreading in Nature, Technology and Society. Springer, 2018. [2] L. Chen, K. Painter, C. Surulescu, and A. Zhigun, “Mathematical models for cell migration: a non-local perspective,” Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 375, p. 20190379, Sept. 2020. 13 (a) (b) FIG. 5.4: Decay of fmax for the diffusion of a particle along a path graph of N = 1001 nodes using the GDE with periodic temporal alternancy, with s1 = 200, s2 = 2 and α = 0.5 (panel (a)) and α = 0.9 (panel (b)). The solid (red) curve is the best power-law fitting for the fmax vs. t for 0 < t < 2, and the broken (blue) line is the same for 2 < t < 4 (see text for details). [3] A. E. B. T. King and M. S. Turner, “Non-local interactions in collective motion,” Royal Society Open Science, vol. 8, pp. rsos.201536, 201536, Mar. 2021. [4] C. Yu, J. Guan, K. Chen, S. C. Bae, and S. Granick, “Single-Molecule Observation of Long Jumps in Polymer Adsorption,” ACS Nano, vol. 7, pp. 9735–9742, Nov. 2013. [5] R. Metzler and J. Klafter, “The random walk’s guide to anomalous diffusion: a fractional dynamics approach,” Physics Reports, vol. 339, pp. 1–77, Dec. 2000. [6] I. M. Sokolov, “Models of anomalous diffusion in crowded environments,” Soft Matter, vol. 8, no. 35, pp. 9043–9052, 2012. [7] B. Li and J. Wang, “Anomalous heat conduction and anomalous diffusion in one-dimensional systems,” Physical review letters, vol. 91, no. 4, p. 044301, 2003. [8] M. Kong and B. Van Houten, “Rad4 recognition-at-a-distance: Physical basis of conformation-specific anomalous diffusion of dna repair proteins,” Progress in biophysics and molecular biology, vol. 127, pp. 93–104, 2017. [9] M. Barbi, C. Place, V. Popkov, and M. Salerno, “A model of sequence-dependent protein diffusion along dna,” Journal of biological physics, vol. 30, no. 3, pp. 203–226, 2004. [10] L. Liu, A. G. Cherstvy, and R. Metzler, “Facilitated diffusion of transcription factor proteins with anomalous bulk diffusion,” The Journal of Physical Chemistry B, vol. 121, no. 6, pp. 1284–1289, 2017. [11] M. Weiss, M. Elsner, F. Kartberg, and T. Nilsson, “Anomalous subdiffusion is a measure for cytoplasmic crowding in living cells,” Biophysical journal, vol. 87, no. 5, pp. 3518–3524, 2004. [12] D. S. Banks and C. Fradin, “Anomalous diffusion of proteins due to molecular crowding,” Biophysical journal, vol. 89, no. 5, pp. 2960–2971, 2005. [13] I. Golding and E. C. Cox, “Physical nature of bacterial cytoplasm,” Physical review letters, vol. 96, no. 9, p. 098102, 2006. [14] S. Gupta, R. Biehl, C. Sill, J. Allgaier, M. Sharp, M. Ohl, and D. Richter, “Protein entrapment in polymeric mesh: Diffusion in crowded environment with fast process on short scales,” Macromolecules, vol. 49, no. 5, pp. 1941–1949, 2016. [15] P. Tan, Y. Liang, Q. Xu, E. Mamontov, J. Li, X. Xing, and L. Hong, “Gradual crossover from subdiffusion to normal diffusion: a many-body effect in protein surface water,” Physical review letters, vol. 120, no. 24, p. 248101, 2018. [16] C. Jiang, C. Cui, L. Li, and Y. Shao, “The anomalous diffusion of a tumor invading with different surrounding tissues,” PloS one, vol. 9, no. 10, p. e109784, 2014. [17] P. Bursac, G. Lenormand, B. Fabry, M. Oliver, D. A. Weitz, V. Viasnoff, J. P. Butler, and J. J. Fredberg, “Cytoskeletal remodelling and slow dynamics in the living cell,” Nature materials, vol. 4, no. 7, pp. 557–561, 2005. [18] N. Shimamoto, “One-dimensional diffusion of proteins along dna: its biological and chemical significance revealed by single-molecule measurements,” Journal of Biological Chemistry, vol. 274, no. 22, pp. 15293–15296, 1999. 14 FIG. 5.5: Decay of fmax for the diffusion of a particle along a path graph of N = 1001 nodes using the GDE with periodic temporal alternancy, with s1 = 200, s2 = 2 and α = 0.5 (blue circles) and α = 0.9 (red triangles). The solid lines are used to guide the eye. The x-axis is in logarithmic scale to maximize the visualization effects. [19] J. Gorman and E. C. Greene, “Visualizing one-dimensional diffusion of proteins along dna,” Nature structural & molecular biology, vol. 15, no. 8, pp. 768–774, 2008. [20] J.-J. Song, R. Bhattacharya, H. Kim, J. Chang, T.-Y. Tang, H. Guo, S. K. Ghosh, Y. Yang, Z. Jiang, H. Kim, et al., “One-dimensional anomalous diffusion of gold nanoparticles in a polymer melt,” Physical review letters, vol. 122, no. 10, p. 107802, 2019. [21] Y. Sagi, M. Brook, I. Almog, and N. Davidson, “Observation of anomalous diffusion and fractional self-similarity in one dimension,” Physical review letters, vol. 108, no. 9, p. 093002, 2012. [22] D. Stauffer, C. Schulze, and D. W. Heermann, “Superdiffusion in a model for diffusion in a molecularly crowded environment,” Journal of biological physics, vol. 33, no. 4, p. 305, 2007. [23] G. I. Livshits, A. Stern, D. Rotem, N. Borovok, G. Eidelshtein, A. Migliore, E. Penzo, S. J. Wind, R. Di Felice, S. S. Skourtis, et al., “Long-range charge transport in single g-quadruplex dna molecules,” Nature nanotechnology, vol. 9, no. 12, pp. 1040–1046, 2014. [24] H. G. Schmidt, S. Sewitz, S. S. Andrews, and K. Lipkow, “An integrated model of transcription factor diffusion shows the importance of intersegmental transfer and quaternary protein structure for target site finding,” PLOS one, vol. 9, no. 10, p. e108575, 2014. [25] M. Sheinman and Y. Kafri, “The effects of intersegmental transfers on target location by proteins,” Physical biology, vol. 6, no. 1, p. 016003, 2009. [26] D. Krepel and Y. Levy, “Intersegmental transfer of proteins between dna regions in the presence of crowding,” Physical Chemistry Chemical Physics, vol. 19, no. 45, pp. 30562–30569, 2017. [27] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. Hwang, “Complex networks: Structure and dynamics,” Physics Reports, vol. 424, pp. 175–308, Feb. 2006. [28] E. Estrada, The Structure of Complex Networks: Theory and Applications. Oxford: Oxford University Press, 2011. [29] V. E. Tarasov, “No nonlocality. no fractional derivative,” Communications in Nonlinear Science and Numerical Simulation, vol. 62, pp. 157–163, 2018. [30] M. Du, Z. Wang, and H. Hu, “Measuring memory with the order of fractional derivative,” Scientific reports, vol. 3, no. 1, pp. 1–3, 2013. [31] E. Estrada, “Path Laplacian matrices: Introduction and application to the analysis of consensus in networks,” Linear Algebra and its Applications, vol. 436, pp. 3373–3391, May 2012. 15 [32] E. Estrada, E. Hameed, N. Hatano, and M. Langer, “Path Laplacian operators and superdiffusive processes on graphs. I. One-dimensional case,” Linear Algebra and its Applications, vol. 523, pp. 307–334, June 2017. arXiv: 1604.00555. [33] E. Estrada, E. Hameed, M. Langer, and A. Puchalska, “Path Laplacian operators and superdiffusive processes on graphs. II. Two-dimensional lattice,” Linear Algebra and its Applications, vol. 555, pp. 373–397, Oct. 2018. [34] E. Estrada, “Path Laplacians versus fractional Laplacians as nonlocal operators on networks,” New Journal of Physics, vol. 23, p. 073049, July 2021. [35] V. Balakrishnan, “Anomalous diffusion in one dimension,” Physica A: Statistical Mechanics and its applications, vol. 132, no. 2-3, pp. 569–580, 1985. [36] K. Anguige and C. Schmeiser, “A one-dimensional model of cell diffusion and aggregation, incorporating volume filling and cell-to-cell adhesion,” Journal of mathematical biology, vol. 58, no. 3, pp. 395–427, 2009. [37] D. Villamaina, A. Sarracino, G. Gradenigo, A. Puglisi, and A. Vulpiani, “On anomalous diffusion and the out of equilibrium response function in one-dimensional models,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2011, no. 01, p. L01002, 2011. [38] J. Padgett, E. Kostadinova, C. Liaw, K. Busse, L. Matthews, and T. Hyde, “Anomalous diffusion in one-dimensional disordered systems: a discrete fractional laplacian method,” Journal of Physics A: Mathematical and Theoretical, vol. 53, no. 13, p. 135205, 2020. [39] S. Nakade, K. Kanki, S. Tanaka, and T. Petrosky, “Anomalous diffusion of a quantum brownian particle in a onedimensional molecular chain,” Physical Review E, vol. 102, no. 3, p. 032137, 2020. [40] S. Jespersen, R. Metzler, and H. C. Fogedby, “Levy flights in external force fields: Langevin and fractional Fokker-Planck equations and their solutions,” Physical Review E, vol. 59, pp. 2736–2745, Mar. 1999. [41] B. Dybiec, E. Gudowska-Nowak, E. Barkai, and A. A. Dubkov, “Levy flights versus Levy walks in bounded domains,” Physical Review E, vol. 95, p. 052102, May 2017. [42] J. Klafter, A. Blumen, and M. F. Shlesinger, “Stochastic pathway to anomalous diffusion,” Physical Review A, vol. 35, pp. 3081–3085, Apr. 1987. [43] A. Allen-Perkins, A. B. Serrano, T. A. de Assis, and R. F. S. Andrade, “Approach to the inverse problem of superdiffusion on finite systems based on time-dependent long-range navigation,” Physical Review E, vol. 100, p. 030101, Sept. 2019. [44] R. Murugan, “Generalized theory of site-specific dna-protein interactions,” Physical Review E, vol. 76, no. 1, p. 011901, 2007. [45] E. F. Koslover, M. D. de la Rosa, and A. J. Spakowitz, “Crowding and hopping in a protein’s diffusive transport on dna,” Journal of Physics A: Mathematical and Theoretical, vol. 50, no. 7, p. 074005, 2017. [46] A. Reynolds, “On the anomalous diffusion characteristics of membrane-bound proteins,” Physics Letters A, vol. 342, no. 5-6, pp. 439–442, 2005. [47] R. Garrappa and M. Popolizio, “Computing the Matrix Mittag-Leffler Function with Applications to Fractional Calculus,” Journal of Scientific Computing, vol. 77, pp. 129–153, Oct. 2018. [48] S. B. Alves, G. F. de Oliveira, L. C. de Oliveira, T. Passerat de Silans, M. Chevrollier, M. Oria, and H. L. de S. Cavalcante, “Characterization of diffusion processes: Normal and anomalous regimes,” Physica A: Statistical Mechanics and its Applications, vol. 447, pp. 392–401, Apr. 2016.