A space-fractional reaction-diffusion system with cylindrical symmetry

A space-fractional reaction-diffusion system with cylindrical symmetry Dimiter Prodanov ∗ ∗Lab. of Neurotechnology, PAML-LN, IICT, Bulgarian Academy of Sciences, 1431 Sofia, Bulgaria (e-mail: dimiter.pro[email protected]as.bg). Abstract: Diffusion within porous media, such as biological tissues, exhibits departures from conventional Fick’s laws, which could result in space-fractional diffusion. The paper considers a reaction-diffusion system with two spatial compartments – a proximal one of finite radius having a source, and an outer one extending to infinity where the source is not present but firstorder decay of the diffusing species takes place. The system models the foreign body reaction around an implanted electrode. Microscopic heterogeneity inside the tissue was modeled by a space-fractional Riesz Laplacian acting on the concentration. This allows for a flexible approach when estimating transport parameters from experimental data. The steady-state of the system is solved in terms of Hankel and Mellin transforms, resulting in a Fox H-function. In the integerorder case, the analytical solution reduces to a superposition of modified Bessel functions of the first and second kinds. Solutions are exhibited by numerical quadrature of the involved Bessel function integrals. Keywords: Hankel transform, Mellin transform, Riesz Laplacian, Fox H-function, Bessel function 1. INTRODUCTION The diffusion in porous media, such as biological tissues, is characterized by deviations from the usual Fick’s diffusion laws (see discussion in Postnikov et al. (2022); Metzler (2020)). One such type of anomalous diffusion is the spacefractional one that arises as an appropriate limit of the Continuous Time Random Walk (CTRW), when jump lengths follow a L´evy stable distribution (Metzler and Klafter, 2000). The source-free diffusion equation can be derived from a non-local generalization of the flux j=D∇βc where Dis the generalized diffusion constant, and c is the concentration of the diffusing species, by using a co-ordinate independent fractional gradient operator ∇β(ˇ Silhav´y, 2019). Microscopic heterogeneity inside the tissue can be lumped into the fractional order parameter β. This potentially allows for a more flexible approach when estimating transport parameters from experimental data. On unbounded domain, mass conservation results in the continuity equation ∂tc=∇·j=−D(−∆)αc, α =1 + β 2(1) where the symbol (−∆)αdenotes the Riesz fractional Laplacian operator of order α(Riesz, 1938). This formulation informs the physical interpretation of the spacefractional diffusion equation. 2. THE REACTION-DIFFUSION SYSTEM A one-dimensional reaction-diffusion application, modeling the spatial distribution of reacting species, was introduced in Prodanov and Delbeke (2016), while the numerical aspects of the two dimensional case have been investigated in Prodanov (2022). The present paper considers a first-order reaction-diffusion system in two dimensions ∂tc=−D(−∆)αc+σ−q c, 0< α ≤1 (2) where σis the intensity of the spatially-extended source and qis an elimination rate. The system can support a nontrivial steady state, which will be studied in the present contribution. Si=σ−qc So=−qc L Fig. 1. Geometry of the reaction-diffusion system We will assume cylindrical geometry of infinite extent, which effectively renders the problem 2 dimensional. The spatial arrangement is represented in Fig. 1. The system is compartmentalized into two spatial compartments – a proximal one of radius L, having a source of intensity σ; and an outer one, which extends to infinity and where the source is not present but there is only a first-order decay. The source term S, therefore splits into two equations Si=σ−qc in the proximal compartment and So=−qc in the distal one. From now on we denote the steady state concentration with the same label. The steady state of the system (eq. 2) can be written as −D(−∆)αc+σ(r)−qc = 0 (3) with boundary condition c(∞) = 0. For simplicity of the presentation we assume D= 1. If D= 1 then the equation can be reparametrized as σ′=σ/D,q′=q/D. 3. INTEGER-ORDER CASE Since the integer-order Laplacian is local the system can split into two disjoint compartments, which interact only at their mutual boundary. The full solution can be written as c(r) = 1(L−r)cp(r) + 1(r−L)cd(r) (4) where 1() denotes the unit step function under the convention 1(0) = 1/2. The equation for the proximal compartment becomes −(−∆)αcp+σ−qcp= 0 (5) while for the outer one the original system becomes the eigen system −(−∆)αcd−qcd= 0 (6) In the present context, Lis the radius of the source, having dimension of length. The solution for the integer-order case α= 1 can be obtained directly by expressing the Laplacian in cylindrical coordinates. One obtains the ordinary differential equation 1 rc′ p+rc′′ p−q cp=−σ(7) where the prime denotes partial derivation in r. The lefthand side of the equation can be recognized as the modified Bessel equation of order 0. Therefore, for the proximal compartment it holds cp(r) = σ/q +k1I0(√qr) (8) where the coefficient k1is determined from the boundary condition at r=L. The solution for the outer compartment then simply is cd=k2K0(√qr) (9) for the indeterminate coefficient k2. In the above equations, I0and K0denote the modified Bessel functions of order 0. To ensure smoothness over the boundary between compartments at r=Lwe impose matching conditions: cp(L) = cd(L), c′ p(L) = c′ d(L) This results in a linear system for the coefficients, which can be solved as k1=−K1L√qσ Pq ,(10) k2=I1L√qσ Pq (11) where P=K0(L√q)I1(L√q) + I0(L√q)K1(L√q) = 1 L√q(12) using the Wronskian identity. Therefore, the final form of the integer-order solution becomes cp(r) = σ q−σL √qK1(√qL)I0(√qr) (13) cd(r) = σL √qI1(√qL)K0(√qr) (14) The solution is exhibited in Fig. 2. The full solution (eq. 16) is compared to the inner (eq. 13) and outer (eq. 14) solutions for α= 1. Fig. 2. Comparison of the integer-order solutions 4. FRACTIONAL-ORDER CASE For the fractional-order case, the compartmentalized approach is not directly applicable but the solution can be obtained in the Fourier domain. Denoting |k|=ρ, the stationary system can be transformed as −ρ2αˆc+σLJ1(ρL) ρ−qˆc= 0 (15) from where one obtains the solution in transcendentalalgebraic form ˆc=σL J1(ρL) ρ(ρ2α+q) Since we assume axially-symmetric geometry, the solution in the spatial domain can be obtained by a Hankel transform (Piessens, 2010): c(r) = σL Z∞ 0 J0(ρ r)J1(ρL) ρ2α+qdρ (16) From formal perspective this completes the analysis of the system. However, the numerical inversion of the Hankel transform is a demanding problem. Therefore, we will look for analytical approximations of the integral. 4.1 An asymptotic solution in terms of H-functions Remarkably, in the integer-order case the solution for the outer compartment can also be obtained in a different way. Notably, we assume a fictitious delta source of intensity σ′ on the boundary r=L. For such an ring-like source the stationary equation reads ∆cd−qcd+σ′δ(r−L) r= 0 (17) Therefore, in the Hankel domain one obtains −ρ2ˆc−qˆc=−σ′J0(ρL) =⇒ˆc=σ′J0(ρL) ρ2+q Therefore, cd(r) = σ′Z∞ 0 J0(ρ r)J0(ρ L)ρ ρ2+qdρ (18) holds in the spatial domain. In this case σ′/σ =LI1/I2, where Idenotes the value of the respective Bessel integrals at L. The matching can be appreciated in Fig. 2. In the fractional-order case, due to the non-local character of the Riesz Laplacian this approximation does not recover the solution but achieves only an asymptotic similarity. The differences in the obtained solutions can be appreciated in Fig. 3. Eq. 18 can be used as a definition of an asymptotic solution also in the fractional case where in the kernel denominator the exponent 2 is substituted by 2α. The full solution (eq. 16) is compared to the outer one (eq. 18) for α= 8/9. Fig. 3. Comparison of the fractional-order solutions A more tractable asymptotic solution can be obtained when one moves the impulse source to the origin – i.e. L= 0. In this case ˆc=σ |k|2α+q(19) is obtained as a solution in the Fourier domain. Therefore, the asymptotic solution is obtained as the Hankel transform: ca(r) = σZ∞ 0 J0(ρ r)ρ ρa+qdρ, a = 2α(20) The function is plotted in Fig. 6 with αparametrization. Remarkably, for a= 2 for the outer component we obtain c2(r) = σ′Z∞ 0 J0(ρ r)ρ ρ2+qdρ =σ′K0(√qr) (21) where, K0is the corresponding modified Bessel function and σ′=k2for this case. In such way, the fractional problem leads to a special function, related to the Bessel K0function (Prodanov, 2022). Furthermore, the asymptotic solution can be represented as Fox H-function by virtue of the Mellin transforms theory Fig. 4. Influence of the fractional exponent on the shape of the full solution (Appendix C). The Mellin transform of the asymptotic solution is given by the integral C(s) = Mx[ca](s) = Z∞ 0 xs−1cd(x)dx = Z∞ 0 dρ Z∞ 0 xs−1J0(ρ x)ρ ρa+qdx (22) We use the known Mellin transform pair of the Bessel J0(x) function (see for example Oberhettinger (1974)) Mx[J0(ρx)](s) = 2s−1 ρs Γ (s/2) Γ (1 −s/2) to obtain C(s) = Z∞ 0 Γs 22s−1ρ1−s Γ1−s 2(ρa+q)dρ = Γs 22s−1 Γ1−s 2Z∞ 0 ρ1−s ρa+qdρ, (23) The last integral can be recognized as an Euler Beta integral: Z∞ 0 ρ1−s ρa+qdρ =q(2−s)/a−1 aB1−2−s a,2−s a Therefore, C(s) = q2−s a−1Γ1−2−s aΓ2−s aΓs 22s−1 aΓ1−s 2(24) The above result allows one to invert the expression by the Mellin-Barnes integral and the residue Theorem. The s-independent pre-factor is given by λ=a pq2/a. The equation defines a Fox H -function with kernel Hm,n p,q (s) = Γs 2Γ2 a−s aΓ1−2 a+s a2s−1 Γ1−s 2(25) From where we read off the parameters m= 2, n= 1, p= 1, q= 3. Therefore, up to a constant factor λ, the inverse Mellin transform is given by the H-function ca(r) = λH2,1 1,3    a √qr 2  1−2 a,1 a 0,1 2,1−2 a,1 a0,1 2    (26) Fig. 5. Asymptotic behavior of ca(z) for α= 0.995 where a=β+ 1 and λ= a √q2 aq . For the given kernel one may take evaluation contour Lto be a Bromwich-type vertical line Re(s)=c with −2< c < 2/a for (a > 1), indented to avoid poles, which places the poles of Γ(s/2) and Γ(1 −2/a +s/a) to the right and Γ(2/a −s/a) poles to the left. It should be noted that for rational values of th exponent athe H function is reducible to a Meijer G-function (see Appendix D) and eventually represented by a finite combination of elementary and special functions. This will be a subject of further studies. For example, let a=2. Then two of the gamma factors cancel and one obtains the simpler kernel H= 2s−1Γ2s 2(27) which is the Mellin transform pair of K0(z), thus confirming eq. 21. In G-function notation K0(z) = 1 2G2,0 0,2z2 4− 0,0(28) Finally, so-identified H-function from eq. 26 allows one to derive asymptotics for large values reven without computing explicitly the function. For the range of interest a∈(1,2] cahas a purely algebraic large-r expansion generated by the right–half-plane poles located at 2 − s=−ak, k ∈N. Therefore, the leading residue (i.e. for k= 1) is c1= Γ a 2+ 122a+1 sin πa 2 πza+2 , a = 2α(29) A plot for α= 0.995 is presented in Fig. 5. 5. NUMERICAL EXPERIMENTS Numerical experiments were performed in the computer algebra system Maxima v. 5.47.0 running in a Windows (TM) 10 operation system. To this end, the DoubleExponential (DE) routine was ported to Lisp and integrated into the computer algebra system Maxima. The code can be downloaded from https://github.com/ dprodanov/intde. Plots of the solutions have been computed and rendered using the DE integration routine. Everywhere q= 1 was used. Fig. 6. Comparison of ca(z) for α=0.85, 0.995 with K0(z) 5.1 The Double-Exponential quadrature method The DE quadrature integration (Takahasi and Mori, 1974; Mori, 1985) can be summarized as follows. The integral on the interval A I=ZA f(x)dx =ZA f(x)dx =Z∞ −∞ f◦ϕ(t)ϕ′(t)dt (30) is transformed to an integral on the entire real line and the function ϕguarantees double exponential convergence to Ias the limit of the Riemannian sum I= lim N→∞ SN= lim N→∞ h N X k=−N ϕ′(kh) | {z } wk f(ϕ(kh) |{z} xk ) (31) Therefore, the value of the integral can be approximated by the truncated sum SN≈I, where the weights are given by wkand the abscissas (i.e. evaluation points) by xk, while his an adjustable parameter. For the case of a semifinite interval the DE method employs the transformation x=ϕ(t) = exp π 2sinh t(32) The absolute error of the method is of the order exp (–c N/ log N). Oscillatory kernels of the form f(x)∼ g(x) sin(ωx +θ) at infinity are supported by the modification of DE method (Ooura and Mori, 1991) using x=ϕ(t) = t 1−exp (−6 sinh t)(33) The outer asymptotic solution ca(r) (eq. 20) is compared for two values of α– 0.85 and 0.995 – with the integerorder one in Fig. 6. 5.2 Convergence acceleration method for numerical Hankel transforms A way to overcome the slow convergence of the Bessel integrals is to evaluate intermediate integrals between the zeros of the Bessel J function (Lucas and Stone, 1995) and to apply convergence acceleration techniques on the partial sums. Even more economic approach is to use the approximation of the Bessel J zeros which can be computed as rν(k) = πk+ν 2−1 4+O(1/k) (34) for large arguments; here νdenotes the order of the Bessel function. An even improved estimate from the same reference is rν(k) = b−(4n2−1)/(8b) + O(1/k2), b=πk+ν 2−1 4(35) Then the Hankel transform of Fcan be evaluated for some large integer Nas the sum f(z) = N X k=0 Zrν(k+1) rν(k) F(r)Jν(rz)rdr+ Z∞ rν(N+1) F(r)Jν(rz)rdr (36) where rν(0) = 0 is formally adjoined to the zero set. The last term of the sum can be used as an error estimate or, alternatively, it can be approximated using the Bessel function asymptotic. A method based on the DE transformation and employing Wynn’s ϵ-algorithm for convergence acceleration has been introduced in Prodanov (2022). An acceptable approximation was found to use only the first 15 asymptotic zeros. The algorithm was employed to plot Fig. 6. 5.3 Plots of the solutions Eqs. 16 and 18, representing the solutions, can be computed numerically to a predefined order of precision. Plots of the solutions are presented in Figs. 2– 4 for L= 1, s= 1, and q= 1 using the DE method for oscillatory integrals. The full solution is plotted for α= 1,3/4,2/3 in Fig. 4. The impact of the fractional exponent can be appreciated in Figs. 5 and 6. From the plot the heavy trails of the concentration can be appreciated against the baseline of K0(i.e. integer-order solution). 6. DISCUSSION The fractional reaction-diffusion system (eq. 2) has been previously used to model the distribution of diffusing species around an implanted electrode (Prodanov and Delbeke, 2016). Such electrodes are routinely used in neurophysiological studies and deep brain stimulation for sensing neural activity or deep brain stimulation for therapeutic applications. Typical implanted electrodes have circular or rectangular cross-section shapes and high aspect ratios. The fractional Laplacian can be defined on unbounded domains in several equivalent ways (Kwa´snicki, 2017). However, when these definitions are restricted to bounded domains, the associated boundary conditions lead to distinct operators (discussion in Lischke et al. (2020)). Although reasoned through the spectral representation (Appendix B), the present paper favors the use of a potential approach implementing long-range spatial interactions since it connects with the physical interpretation of the diffusion equation. Although the presented paper employs the Riesz Laplacian, the modeling approach is not restricted by the choice of operator and can be extended to other definitions. The present contribution obtained an analytical solution of the problem, which is computable by quadrature (eq. 16) as demonstrated in Fig. 4. From theoretical perspective, of some interest is also the asymptotic solution (eq. 26) obtained in terms of a Fox H-function. Computation of such functions is an open area of research since they are very general objects. Finally, of practical interest is also eq. 29 (see also Fig. 5). Its application would allow for a more flexible approach when estimating transport parameters from experimental data. ACKNOWLEDGEMENTS The present work is funded by the Horizon Europe project VIBraTE, grant agreement no. 101086815. Appendix A. INTEGRAL TRANSFORMS A.1 Fourier transform The Fourier transform will be defined under the ”engineering” convention ˆ f(k) = Ff(x) := ZRd f(x)e−ik·xdxd with an inverse f(x) = F−1ˆ f(k) := 1 (2π)dZRd ˆ f(k)eik·xdxd A.2 Hankel transform convention The automorphic Hankel transform is defined as (Piessens, 2010): ˆ f(ρ) = Hνf(z) := Z∞ 0 f(z)Jν(ρ z)zdz (A.1) Remarkably, for 2 spatial dimensions the Laplacian is represented by a monomial factor: H0∆f(z) = −ρ2ˆ f(ρ) and in the similar way the image of the Riesz Laplacian is (−∆)α7→ −|ρ|2α. A.3 Mellin transform The Mellin transform of a function f(t) is defined as Mt[f](s) := Z∞ 0 ts−1f(t)dt =F(s),(A.2) wherever the integral exists on the complex plane. The inverse Mellin transform is given by complex integration along a Bromwitch contour: f(z) = M−1[F](z) = 1 2πi ZBr F(s)z−sds Appendix B. THE RIESZ LAPLACIAN OPERATOR The Riesz operator can be defined in the Fourier domain (Kwa´snicki, 2017) where −(−∆)αf(x)7→ |k|2αˆ f(k). Then starting from the algebraic substitution (−∆)α7→ −|k|2α, the fractional Laplacian can be interpreted as the gradient of another operator, that is −|k|2α=ik·ik0|k|2α−1,k0=k k where the dot denotes the scalar product and k0is a unit wave vector in the Fourier space. Therefore, the associated Riesz gradient can be defined as a pseudo-differential operator ∇β7→ ik0|k|1−β=ik/|k|β, β = 2α−1 for a suitable function space. This has the advantageous physical interpretation of a generalized first Fick’s law, which This corresponds to the convolution with ∇β=I1−β∗∇ where I1−βdenotes the Riesz potential. Appendix C. FOX H FUNCTIONS The Fox H-function Hm,n p,q [z|·] is defined by a Mellin-Barns integral as follows: Hm,n p,q z (a1, A1),(a2, A2),...,(am, Am) (b1, B1),(b2, B2),...,(bn, Bn):= 1 2πi ZLHm,n p,q (s)·z−sds =1 2πi ZLQm j=1 Γ(bj+Bjs) Qq j=1+mΓ(1 −bj−Bjs)· Qn j=1 Γ(1 −aj−Ajs) Qp j=n+1 Γ(aj+Ajs)·z−sds, (C.1) where: •mand nare non-negative integers, •pand qare positive integers, •zis a complex variable, •aj, bj, Aj, Bjare parameters with positive real parts, whereas Lis a suitable contour in the complex plane that separates the poles of the Gamma functions in the numerator. Appendix D. MEIJER G FUNCTIONS The Meijer G functions are closely related to the Fox Hfunctions. The G function Gm,n p,q [z|·] is defined by a MellinBarns integral as follows: Gm,n p,q z a1, a2, . . . , ap b1, b2, . . . , bq:= 1 2πi ZLGm,n p,q (s)·z−sds =1 2πi ZLQm j=1 Γ(bj−s) Qq j=1+mΓ(1 −bj+s)· Qn j=1 Γ(1 −aj+s) Qp j=n+1 Γ(aj−s)·z−sds, where: •mand nare non-negative integers, •pand qare positive integers, •zis a complex variable, •aj, bjare parameters with positive real parts, and and the contour Lseparates poles. REFERENCES Kwa´snicki, M. (2017). Ten equivalent definitions of the fractional laplace operator. Fract. Calc. Appl. Anal., 20(1). doi:10.1515/fca-2017-0002. Lischke, A., Pang, G., Gulian, M., Song, F., Glusa, C., Zheng, X., Mao, Z., Cai, W., Meerschaert, M.M., Ainsworth, M., and Karniadakis, G.E. (2020). What is the fractional Laplacian? A comparative review with new results. J. Comput. Phys., 404, 109009. doi:10. 1016/j.jcp.2019.109009. Lucas, S. and Stone, H. (1995). Evaluating infinite integrals involving bessel functions of arbitrary order. J. Comput. Appl. Math., 64(3), 217–231. doi:10.1016/ 0377-0427(95)00142-5. Metzler, R. (2020). Superstatistics and non-gaussian diffusion. Eur. Phys. J. Special Topics, 229(5), 711–728. doi:10.1140/epjst/e2020-900210-x. Metzler, R. and Klafter, J. (2000). The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Physics Reports, 339(1), 1–77. doi:10.1016/ s0370-1573(00)00070-3. Mori, M. (1985). Quadrature formulas obtained by variable transformation and the DE-rule. Journal of Computational and Applied Mathematics, 12-13, 119–130. doi:10.1016/0377-0427(85)90011-1. Oberhettinger, F. (1974). Tables of Mellin Transforms. Springer Berlin / Heidelberg, Berlin, Heidelberg. Ooura, T. and Mori, M. (1991). The double exponential formula for oscillatory functions over the half infinite interval. J. Comput. Appl. Math., 38(1?3), 353–360. doi: 10.1016/0377-0427(91)90181-i. Piessens, R. (2010). Hankel Transform. In A. Poularikas (ed.), Transforms and Applications Handbook, Electrical Engineering Handbook ; 43. CRC Press, Boca Raton, Fla, 3rd ed edition. Postnikov, E.B., Lavrova, A.I., and Postnov, D.E. (2022). Transport in the brain extracellular space: Diffusion, but which kind? International Journal of Molecular Sciences, 23(20), 12401. doi:10.3390/ijms232012401. Prodanov, D. (2022). First-Order Reaction-Diffusion System with Space-Fractional Diffusion in an Unbounded Medium, 65–70. Springer International Publishing. doi: 10.1007/978-3-030-97549-4 7. Prodanov, D. and Delbeke, J. (2016). A model of spacefractional-order diffusion in the glial scar. J. Theor. Biology, 403, 97–109. doi:10.1016/j.jtbi.2016.04.031. Riesz, M. (1938). Int´egrales de Riemann–Liouville et potentiels. Acta Sci. Math. Szeged, 9, 1–42. ˇ Silhav´y, M. (2019). Fractional vector analysis based on invariance requirements (critique of coordinate approaches). Continuum Mech. Thermodyn., 32(1), 207– 228. doi:10.1007/s00161-019-00797-9. Takahasi, H. and Mori, M. (1974). Double exponential formulas for numerical integration. Pub. RIMS Kyoto Univ, 9, 721–741.