Reaction-subdiffusion front propagation in a comblike model of spiny dendrites

Fractional reaction-diffusion equations are derived by exploiting the geometrical similarities between a comb structure and a spiny dendrite. In the framework of the obtained equations, two scenarios of reaction transport in spiny dendrites are explored, where both a linear reaction in spines and nonlinear Fisher-Kolmogorov-PetrovskiiPiskunov reactions along dendrites are considered. In the framework of fractional subdiffusive comb model, we develop a Hamilton-Jacobi approach to estimate the overall velocity of the reaction front propagation. One of the main effects observed is the failure of the front propagation for both scenarios due to either the reaction inside the spines or the interaction of the reaction with the spines. In the first case the spines are the source of reactions, while in the latter case, the spines are a source of a damping mechanism.


I. INTRODUCTION
A general scheme of front propagation in reaction-transport equations is attracting much attention that is related to the considerable progress in our understanding of fractional reactionsubdiffusion systems [1][2][3][4][5][6][7][8][9].New experimental findings in calcium transport and reaction transport in neuroscience [10] (see also [11]) open new questions on the geometric impact on the calcium transport in spiny dendrites and extension of a general scheme of this description of reaction under subdiffusion.An important new trend here is investigation of the influence of dendrite geometry on the reaction transport characteristics.
In this paper we consider front propagation in spiny dendrites as the reaction-subdiffusion scheme, where spines play a crucial role in both the reaction and subdiffusion.Spiny dendrites are dendrites with lateral small protrusions, called dendritic spines, located on the surface of the dendrite.They can be found on the dendrites of most principal neurons in the brain, including the pyramidal neurons of the neocortex, the medium spiny neurons of the striatum, and the Purkinje cells of the cerebellum.In Fig. 1 we show a schematic picture of a spiny dendrite.A dendritic spine (or spine) is a small membranous protrusion from a neuron's dendrite that typically receives input from a single synapse of an axon.Dendritic spines serve as a storage site for synaptic strength and help transmit electrical signals to the neuron's cell body.Most spines have a bulbous head (spine head size ∼1 μm) and a thin neck (size ∼0.1 μm) that connects the head of the spine to the shaft of the dendrite.The dendrites of a single neuron can contain hundreds or thousands of spines.Their heads have an active membrane, and as a consequence, they can sustain the propagation of an action potential with a rate that depends on the spatial density of spines [12].
Recent experiments together with numerical simulations investigate the calcium transport inside spiny dendrites to understand the role of calcium in signal transmission and neural plasticity [13][14][15][16].Following experimental findings, some theoretical approaches have been suggested with a different direction of exploration of the transport properties of spiny dendrites.A seminal theoretical result on the front propagation of a solitary wave in dendritic spines is suggested in Ref. [17], where the spines are explored as active clusters, while the dendrite is passive.One finding here is the front propagation failure in terms of the separation between regularly separated clusters.The front propagation failure means the overall velocity of the moving reaction front vanishes and a static front solution exists.There are many mechanisms that produce propagation failure.For example, bistable kinetics, spatial heterogeneities, etc. [18,19].
Experimental observations of CaMKII (Ca 2+ -calmodulindependent protein kinase II, a key regulator of the synaptic function) translocation waves [10] opened an important question on theoretical understanding of the reaction front propagation.To study the translocation waves of CaMKII in spiny dendrites, a system of coupled reaction-diffusion equations was proposed [11], where reactions were considered inside dendrites between two chemicals and a Markovian switching with particles concentrated inside spines [11,17,20].The failure of the propagation in terms of the switching rate has been observed numerically [11] as well.
In the present work we suggest an analytical treatment of reaction-transport scenarios in spiny dendrites, where we explore both a linear reaction in spines and a nonlinear reaction along dendrites.The latter is a realization of the Fisher-Kolmogorov-Petrovskii-Piskunov (FKPP) scheme [21,22].As shown in Refs.[10,11], the process of activation is related to the first-order reaction scheme P p + P a → 2P a that is considered in the framework of the FKPP interaction of the order of P p P a , where P p and P a are concentrations of prime and activated CaMKII.We explore the same scheme in the framework of only a single nonlinear reaction transport equation for the activated CaMKII contaminant P ≡ P a with the interaction term in the form of the logistic or FKPP-like reaction C(P ) = CP (1 − P ), which is widely used in reaction-transport equations [1][2][3][4][5][6][7][8][9].To account for the geometric impact on the reaction front propagation, we suggest a fractional subdiffusive comb model [23] amended by the FKPP reaction.Then a hyperbolic scaling of the dynamical variables is applied to arrive at a Hamilton-Jacobi equation and to estimate the overall velocity of the reaction front propagation.One of the main effects observed is the failure of the front propagation in both scenarios due to either the reaction inside spines or the interaction of the reaction with spines.The Hamilton-Jacobi approach is a technically elegant implement for studying the propagation front velocity under general transport mechanisms and is especially useful when propagation takes place in heterogeneous media [24].For the reaction transport equations, it relates to developing an appropriate technique of treating front propagation, where an appropriate hyperbolic scaling of the reaction-transport equation makes it possible to estimate the overall rate of the spreading reaction wave without resolving its shape [25,26].The method of hyperbolic scaling is based on the introduction of a small parameter ε → 0, rescaling of coordinates and time (x,t) → (x/ε,t/ε), and the contaminant's density distribution function (for example, CaMKII).In this case the problem of the wave propagation reduces to the dynamics of the leading edge or the reaction front.Therefore, one analyzes the reaction-transport behavior in the leading edge, where, in the long-range and long-time limits, the detailed shape of the traveling wave is not important when the front invades the unstable state.

II. GEOMETRIC IMPACT ON ANOMALOUS TRANSPORT IN DENDRITES
As shown in previous studies, the geometric nature of spiny dendrites plays an essential role in kinetics [27][28][29][30][31][32].The real distribution of spines along the dendrite, their size, and their shapes are completely random [33], and inside spines, not only the spine necks but also the spine itself act as a transport barrier [28,32,34].Therefore, a reasonable assumption is to consider anomalous diffusion along both the spines and dendrite.So we propose models based on a comblike structure that mimic a spiny dendrite, where the backbone is the dendrite and the fingers (lateral branches) are the spines; see Fig. 2. In this case the dynamics inside fingers corresponds to spines, while the backbone describes diffusion inside dendrites.Note that the comb model is an analog of a one-dimensional (1D) medium where fractional diffusion has been observed and explained in the framework of a so-called continuous time random walk (CTRW) [35][36][37][38][39]. Before discussing the CTRW consideration in the framework of the comb model, let us explain how anomalous diffusion in the comb model relates to the CaMKII transport along the spiny dendrite and how the geometry of the latter relates to the anomalous transport.As noted above the spine cavities behave as traps for the contaminant transport.As follows from a general consideration of a Markov process inside a finite region, the survival probability to find a particle inside a cavity of an arbitrary form with a finite volume decays exponentially with time t (see, for example, [40]), Here τ is the survival time (mean life time), defined by the minimum eigenvalue of the Laplace operator and determined by the geometry of the cavity.For example, in Refs.[31,32], for spines with a head of volume V and a cylindrical spine neck of length L and radius a, the mean life time is τ = LV /π a 2 D, where D is the diffusivity of the spine.Therefore, the mean probability to find a particle inside the spine after time t, averaged over all possible realizations of τ , is the integral where f (τ ) is the distribution function of the survival times τ (recall that the size and shape of spines are random [33]).
In the simplest case, when the distribution is the exponential f (τ ) = (1/τ 0 ) exp(−τ/τ 0 ), one obtains from Eq. ( 1) that the general kinetics is not Markovian.As one easily calculates from Eq. ( 1), the waiting time probability distribution function (pdf) is a stretched exponent, The situation is more interesting when the distribution of the survival times is the power law f (τ ) ∼ 1/τ 1+α (0 < α < 1).In this case the waiting time pdf is the power law as well, ψ(t) ∼ 1/t 1+α , which leads to subdiffusion motion along the dendrite.This result follows from the CTRW theory since all underlying microprocesses are independent Markovian ones with the same distributions [38,39].Now we explain the physical reason for the possible power law distribution ψ(t).At this point we paraphrase some arguments from Ref. [41] (see their Sec.4.2) with the corresponding adaptation to the present analysis.Let us consider the escape from a spine cavity from a potential point of view, where geometrical parameters of the cavity can be related to a potential U .For example, for the simplest case, mentioned above, it is U = V L/πa 2 , which "keeps" a particle inside the cavity, while Dτ 0 plays the role of the kinetic energy, or the "Boltzmann temperature."Therefore, the escape probability from the spine cavity or well is described by the Boltzmann distribution exp(−U/Dτ 0 ).This value is proportional to the inverse waiting or survival time, As mentioned above, the potential U is random and distributed by the Poisson distribution , where U 0 is an averaged geometrical spine characteristic.The probability to find the waiting time in the interval (t,t + dt) is equal to the probability to find the trapping potential in the interval (U,U + dU ), namely, ψ(t)dt = P (U )dU .Therefore, from Eq. ( 2) one obtains ( Here α = Dτ 0 U 0 ∈ (0,1) establishes a relation between the geometry of the dendrite spines and the subdiffusion observed in Refs.[27,28] and supports application of the comb model, which is a convenient implement for analytical exploration of anomalous transport in spiny dendrites in the framework of the CTRW consideration.

III. FRACTIONAL DIFFUSION ON A COMB WITH REACTION
Geometrically, spiny dendrites in three-dimensional (3D) space are completely described by a comb structure in the two dimensions, where the spine density on the cylinder surface is projected on the 1D axis (say the x axis): ρ(x,r = const,φ) → ρ(x).Here ρ(x,r = const,φ) is the spine density, and ρ(x) is the density of the comb fingers.In what follows, we consider ρ(x) = g = const, which is, probably, the most realistic case.Fractional diffusion inside the spines is described by fractional diffusion inside the fingers.Therefore, one considers a twosided comb model, and the starting mathematical point of the phenomenological consideration is the Fokker-Planck equation obtained in [23].Here we present an alternative way for inferring the fractional Fokker-Planck equation on a comb (9).
The comb model is also known as a toy model for a porous medium used for exploration of low-dimensional percolation clusters [36].It should be admitted that the transport exponent, obtained in the previous section as the geometric impact, is α = 1/2 for the CTRW consideration in the framework of the comb model.Usually, anomalous diffusion on the comb is described by the two-dimensional (2D) distribution function P = P (x,y,t), and a special behavior is that the displacement in the x direction is possible only along the structure axis (x axis at y = 0).Therefore, diffusion in the x direction is highly inhomogeneous.Namely, the diffusion coefficient is D xx = Dδ(y), while the diffusion coefficient in the y direction (along fingers) is a constant D yy = D. Therefore, this inhomogeneous diffusion is described by the Fokker-Planck equation in the dimensionless time and coordinates, It is obtained by rescaling with relevant combinations of the comb parameters D and D, such that the dimensionless time and coordinates are D 3 t/ D2 → t and Dx/ D → x, Dy/ D → y/ √ g, respectively [42], and parameter g can be considered the constant density of the fingers.
A variety of interactions inside spines leads to correlated noises in dendritic spines [43].The strong correlations of that noise lead to anomalous (subdiffusive) motion inside the spines.Following a phenomenological description by the CTRW, this subdiffusion is controlled by a waiting time pdf ψ(t) decaying according to the power law (3).Therefore, normal diffusion of the contaminant density P (x,y,t), for example, activated CaMKII, in spines or fingers is replaced by the anomalous transition term where K(t) is the time memory kernel determined by the waiting time pdf in the Laplace domain For subdiffusion, 0 < γ < 1 and Eq. ( 6) yields [5] One should recognize that Eq. ( 5) is a formal expression for the anomalous transport with a very complicated form in the time domain, which, in turn, is very inconvenient for the analytical treatment.Therefore, the comb model can be presented in the Laplace domain.Substituting Eq. ( 5) in Eq. ( 4), then performing the Laplace transform and taking into account Eq. ( 6), one obtains the comb model in the Laplace domain P (s) = L[P (t)], Here P 0 = P (x,y,t = 0) is the initial condition.As noted, the kernel K is problematic for the Laplace inversion since it leads to the appearance of the initial condition.To overcome this obstacle, one multiplies Eq. ( 7) by s α−1 and then perform the Laplace inversion, which yields Amending this equation by the reaction term, one arrives at the integro-differential equation: In what follows we use convenient notations of fractional integro-differentiation (see, for example, [44,45]).Namely, the α-order fractional integral, for 0 < α < 1, is where (α) is a gamma function.Fractional differentiation used here in the Caputo form, is defined by means of the Laplace transform L[f (t)] = f (s) for the fractional integration L[I α t f (t)] = s −α f (s), which yields Owing to this notation, from Eq. ( 9) the equation for P (x,y,t) reads To give a first and brief insight on the problem of the front propagation, let us consider the linear reaction and γ = 1.In this case, one obtains a "simple" solution [46] for the traveling wave along the x axis (inside dendrites).Introducing the total probability distribution function P 1 (x,t) = dyP (x,y,t), one obtains This yields the coordinates of the front x ∼ t 3/4 that spreads with the decaying velocity v ∼ t −1/4 .This solution illustrates the asymptotic failure of the reaction-transport front propagation due to subdiffusion inside spiny dendrites.

IV. PROTEIN KINASE TRANSLOCATION IN A FKPP COMB
Recently, a mechanism of the translocation wave of CaMKII was suggested [11].As shown, activated CaMKII contaminant travels along dendrites with additional translocation inside spines.The process of activation (the conversion of primed CaMKII to its active state) corresponds to the irreversible reaction that, in the absence of spines, is described by the FKPP equation [11].Therefore, in the framework of the above-suggested scheme of the dispersive subdiffusive comb (10), nonlinear reaction at subdiffusion in dendrites takes place along the x-axis bound, while subdiffusion in fingers describes the translocation inside spines.Therefore, reaction-transport equation (10) now reads Here D describes the diffusivity inside dendrites, while C(P ) = CP (1 − P ) is the nonlinear reaction term.Again, integrating over y to obtain the kinetic equation for the total distribution P 1 (x,t), we have For brevity, we denote P 0 = P (x,y = 0,t).Manipulation with the Laplace transform P1 (x,s) = L[P 1 (x,t)] and substituting P0 (x,s) = s γ 4g P1 (x,s) yield, from Eq. ( 13), the following equation for the Laplace image (see the Appendix): Multiplying this equation by e st and using the identity e st s α f (s) = ∂ ∂t e st s α−1 f (s), we integrate with the corresponding contour to obtain the inverse Laplace transform.This yields The nonlinear term is obtained by the following chain of transformations: Note that a specific property of these transformations is an irreversibility with respect to the Laplace transform since, as is well known, the Laplace transform of the Riemann-Liouville fractional derivative involves the (quasi)initial value terms like P 1 (t = 0) = δ(x) [39,44].
To evaluate the overall velocity of the asymptotic front, let us introduce a small parameter, say ε, at the derivatives with respect to time and space [25,26].To this end we rescale x → x/ε and t → t/ε, and P 1 (x,t) → P ε 1 (x,t) = P 1 ( x ε , t ε ).Therefore, one looks for the asymptotic solution in the form of the Green's approximation, also loosely known as a WKB approximation, The main strategy of implication of this construction is the limit ε → 0; one has exp[− G ε (x,t) ε ] = 0, except for the condition when G ε (x,t) = 0.This equation determines the position of the reaction spreading front; see Eq. (11).Moreover, we consider the limit G(x,t) = lim ε→0 G ε (x,t) as the principal Hamiltonian function [25,26], which makes it possible to apply the Hamiltonian approach to calculate the propagation front velocity.In this case partial derivatives of G(x,t) with respect to time and coordinate have the physical senses of the Hamiltonian and momentum: Now the method of the hyperbolic scaling, explained above, can be applied.Therefore, we have ansatz (16) for the probability distribution function inside dendrites.Inserting expression (16) in Eq. ( 15), one considers fractional integrations in time.Let us start from the last term in Eq. (15), which is the reaction term.We rewrite it in the following convenient form: Then performing the expansion and substituting this in Eq. ( 18), one obtains It should be noted that we neglect differentiation of the upper limit of the integral since this term is of the order of O(ε 1+γ /2 ) ∼ o(ε), which vanishes in the limit ε → 0. The same procedure of expansion is performed for the diffusion term in Eq. ( 15), which yields Differentiating in the limit ε → 0 and taking into account that the Hamiltonian H and the momentum p in Eq. ( 17) are independent of x and t explicitly (which leads to the absence of mixed derivatives), one obtains the Laplace transform of the subdiffusive kernel t − γ 2 .After these procedures in Eqs. ( 19) and ( 20), the kinetic equation ( 15) becomes a kind of Hamilton-Jacobi equation that establishes a relation between the Hamiltonian and the momentum, and the action is G(x,t) = t 0 [p(s) ẋ(s) − H (p(s),x(s))]ds.The rate v at which the front moves is determined at the condition G(x,t) = 0. Together with the Hamilton equations, this yields Note that the first equation in Eq. ( 22) reflects the dispersion condition, while the second one is a result of the asymptotically free particle dynamics, when the action is G(x,t) = px − H t.
Taking into account x = vt, one obtains Eq. ( 22) (see also details of this discussion, e.g., in Refs.[26,47,48]).The combination of these two equations can be replaced by We also have from the front velocity conditions ( 22) ∂ ∂p ln H = 1/p, which, eventually, yields from Eq. ( 23) To proceed, we first note that the limiting case of this result with γ = 0 corresponds to the CaMKII propagation along the dendrite only (i.e., there are no lateral branches or fingers).Therefore, Eq. ( 24) after rescaling D/ √ g → D and C/ √ g → C recovers the FKPP scheme for γ = 0, which yields v = √ DC.
The absence of the failure of the activation front propagation should be noted.This has a simple explanation due to the absence of a reaction "sink" term −hP in Eq. ( 12) by neglecting the possibility of spines to accumulate a large amount of Ca 2+ [13,14], where h is a translocation or accumulation rate [11].Introducing this term in Eq. ( 12), our anticipation is that the hyperbolic scaling for this new equation yields a solution similar to Eq. ( 23) with H = 0, which corresponds to the failure of the front propagation.Moreover, this asymptotic solution for P 1 (x,t) always takes place as one of the possible solutions.
Inserting the sink term in Eq. ( 12), one obtains Repeating the same procedures of the Laplace transform of P and integration over y with the definition P1 = ∞ −∞ P (x,y,s)dy and using the substitute Here we neglect the nonlinear term since, as follows from the above analysis, in the further hyperbolic scaling approximation, this term does not contribute to the Hamilton-Jacobi equation, and we also know how to handle it.Again multiplying this equation by e st and using the same identity e st s α f (s) = ∂ ∂t e st s α−1 f (s) as above, we obtain the inverse Laplace transform.Thus, Eq. ( 26) reads Application of the hyperbolic scaling to the asymptotic solution (16) yields Let us consider the specific case γ = 2/3, which yields For C > 2hg 3 2 there is no failure, and the front asymptotically propagates with a constant velocity.For C 2hg 3 2 the only solution is H = 0 and yields v = 0.So 2hg 3 2 is the minimum reaction rate necessary to sustain propagation along the spiny dendrite due to the presence of translocation.Analogously, (C/2h) 2/3 can also be viewed as the minimum value for the density of spines necessary to have propagation failure.Both results are in agreement with the results obtained from very different models based on the cable model [17].
In a general case, one compares the interplay between the activation CH γ 2 and the translocation −2hg 3 2 H 1−γ terms in Eq. ( 28) in the limit H → 0. For γ ∈ [ 2 3 ,1), the translocation term is dominant and leads to the solution with H = 0 and the failure of the front propagation.When 0 < γ < 2  3 , the activation in dendrites can be dominant.This situation is more complicated, and the activation-translocation front can propagate with an asymptotically finite velocity.

V. ANALYTICAL SOLUTION OF EQ. (12)
Let us consider the linear counterpart of Eq. ( 12) with the linear reaction term C(P ) = CP .We rewrite this equation for the total distribution P 1 (x,t).As follows from the fractional differentiation of Eq. ( 15), this equation reads (see also [23]) with the initial condition P 1 (x,t = 0) = δ(x).After the Fourier transform F[P 1 (x,t)] = P1 (k,t), one obtains the solution in the form of the Mittag-Leffler function, where A(k) = (C − Dk 2 )/2 √ g.At the asymptotic condition, when x ,t 1, we have C Dk 2 , which yields asymptotic behavior of the Mittag-Leffler function as a growing exponent (for the large positive argument) [49], After the Fourier inversion, one obtains which, finally, yields the nonzero and constant overall velocity of the reaction front propagation.Note that for normal diffusion, γ = 0, one arrives at the Fisher velocity v = √ DC/g → √ DC; see the limiting case γ = 0 in Eq. (24).

VI. DISCUSSION
In the paper we show that a comb is a convenient model for analytical exploration of anomalous transport in spiny dendrites.Although understanding of the role of calcium in signal transmission and neural plasticity is mainly due to experimental and numerical studies [11,[13][14][15][16], a simple toy model, like the comb model, makes it possible to suggest and understand a variety of reaction-transport schemes, including anomalous transport, by applying the strong machinery of fractional calculus and hyperbolic scaling for asymptotic methods.
In the present consideration we suggest an analytical description of reaction-transport scenarios in spiny dendrites, where we consider both a linear reaction in spines [see Eqs. ( 9) and (10)] and a nonlinear reaction along dendrites, considered in the framework of the FKPP scheme [21,22].To this end we suggest a fractional subdiffusive comb model, where we apply a Hamilton-Jacobi approach to estimate the overall velocity of the reaction front propagation.We proposed an alternative approach of a recently suggested mechanism of the translocation wave of CaMKII [11], where activated CaMKII contaminant travels along dendrites with additional translocation inside spines, and the process of activation corresponds to the irreversible reaction described by the FKPP equation (15).One of the main effects observed in the framework of the considered model is the failure of the front propagation due to either the reaction inside spines or the interaction of the reaction with spines.In the first case the spines are the source of reactions, while in the latter case the spines are a source of damping; for example, they act as a sink of an activated contaminant (CaMKII).The situation is controlled by three parameters, CaMKII activation C, CaMKII translocation rate h, and the fractional transport exponent γ .The latter reflects the geometrical structure of the transport system: when 0 < γ < 2  3 , the activation in dendrites can be dominant, and the activation-translocation front can propagate with an asymptotically nonzero and constant velocity.For γ = 2/3 we have found a criterion for the emergence of propagation failure or for sustaining the propagation in terms of the reaction rate, the translocation rate, and the spine's density.
It should be noted, in conclusion, that the physical arguments suggested above explain why anomalous transport, namely, subdiffusion, of either CaMKII or neutral particles is possible and support the implementation of the comb model.These arguments are based on the geometry of dendritic spines, which determines an expression for the transport exponent in Eq. ( 3).This situation becomes more sophisticated in the case of the nonlinear FKPP reaction.Indeed, as shown, the power law kernel of the transition probability considered due to the geometric arguments is insensitive to the nonlinear reaction.This consideration differs completely from a mesoscopic non-Markovian approach, developed in Refs.[18,50], where spines-dendrite interaction and an extension including reactions in the spines have been described in the framework of variable residence time.This leads to the essential complication of the transition probability due to the nonlinear reaction term [4,18].An interplay between geometry and the residence time is an interesting feature, which remains for future investigations [51].where for the initial condition we take P (t = 0) = δ(x)δ(y).
For the next step, we look for the solution in the form P (x,y,s) = exp[− s γ /g|y|]f (x,s) .(A4)
which describes the 2D inhomogeneous reaction diffusion in the dispersive medium.Here C[P (x,y,t)] ≡ C(P ) is a reaction kinetic term.In particular, to model reaction kinetics inside dendrites, it can be considered either linear, C(P ) = CP , or logistic, C(P ) = CP (x,y,t)[1 − P (x,y,t)][20].Integration with the power law kernel t −γ ensures anomalous diffusion in both the dendrite and spines.