Emergence of Fractional Kinetics in Spiny Dendrites

Fractional extensions of the cable equation have been proposed in the literature to describe transmembrane potential in spiny dendrites. The anomalous behavior has been related in the literature to the geometrical properties of the system, in particular, the density of spines, by experiments, computer simulations, and in comb-like models.~The same PDE can be related to more than one stochastic process leading to anomalous diffusion behavior. The time-fractional diffusion equation can be associated to a continuous time random walk (CTRW) with power-law waiting time probability or to a special case of the Erd\'ely-Kober fractional diffusion, described by the ggBm. In this work, we show that time fractional generalization of the cable equation arises naturally in the CTRW by considering a superposition of Markovian processes and in a {\it ggBm-like} construction of the random variable.


Introduction
Neurons are the fundamental structural units of the nervous system. These cells are specialized to communicate with each other through electrical and chemical signals, specifically called neural signals. Despite the incredible diversity existing between different neuron types, the basic mechanism to exchange electrical signals is the same as for other excitable cells, and is driven by transmembrane ion currents, generating a variation in the transmembrane voltage V m . The cell membrane is composed by a phospholipid bilayer that isolates the inner part of the cell from the surround, crossed by several macromolecular structures (proteins) that allow ions and other molecules to flow in and out. These biological constituents have been historically described by the use of electrical circuit elements to model the transmembrane potential of the cells. A capacitor element was introduced to mimic the role of the phospholipid bilayer, which keeps a different charge concentration in and out the cell. To model the presence of embedded proteins was considered a resistance, in parallel to the capacitor, to roughly describe the flow of charges in and out.
Dendrites and axons are nonisopotentential parts of the neuron in which the membrane can be geometrically approximated by a cylinder , with the longitudinal axis length greater than the radius. Neglecting radial flow phenomena, the cell membrane is described by a linear density of RC modules, composed by a transmembrane capacitance c m and resistance r m in parallel, connected by an internal (and eventually external) resistance r i associated to the ionic flow parallel to the membrane in the cytoplasm viscous medium inside the cell (or the extracellular fluid outside). The differential equation describing the voltage in this circuit system is a diffusion equation with the addition of an extra term that accounts for the decay of the electric signal in time: where, following [Magin(2006)], the adimensional variables X = x/λ and T = t/τ are introduced. The constants λ = r m /r i and τ = r m c m , the space and time scales of the process, are determined by the values of the membrane resistance and capacitance per unit length of the system. The resulting fundamental solution of the Chauchy problem is a Gaussian suppressed by an exponential decay: the appearance of exponential decay is determined by the term −V m in the equation, which accounts in fact for particles lost in the system through the transmembrane channels. In this system, the mean square displacement (MSD) of diffusing molecules is expected to scale linearly with time. The motion of ions and the subsequent propagation of the action potential along the axons is well described by this model; however, anomalous diffusion behavior in the propagation of the subthreshold potential in dendrites has been measured by several experiments [Nimchinsky(2002), Ionescu(2017), Jacobs(1997)]. Experimental evidences of anomalous diffusion of an inert tracer in spiny branches of Purkinje cells [Santamaria(2006)] suggested that the origin of anomalous diffusion in this system was related to the geometry of the system, and that anomalous diffusion was related to the presence of spines more than to the presence of branches. Furthermore, spine density can change dynamically depending on neuronal activity. The idea of correlation between spines density and the anomalous time scaling exponent of the MSD was suggested in [Henry(2008)], because spine density is an important feature for the physiological behavior of several types of neurons, and the low density of spines is associated to aging [Jacobs(1997), Duan(2003], neurological disorders [Nimchinsky(2002)] and syndromes [Suetsugu(1980)].
In [Henry(2008)], the anomalous diffusion of ions was introduced modifying the Nernst-Planck equation (NPE) to generate a fractional Brownian motion (fBm) and a continuous time random walk (CTRW) process. In these cases the diffusion coefficient D is not constant as in the standard NPE, but it becomes a time dependent operator characterized by the scaling parameter 0 < α ≤ 1. For fBm and CTRW respectively we have: where ∂ 1−α ∂t 1−α is the Riemann-Liouville fractional derivative operator. This approach leads to the following differential equations for the two models: and where the terms −µ 2 κT κ−1 V f Bm and −µ 2 ∂ κ−1 VCT RW ∂T κ−1 still account for particles loss in the system instead of the term −V m (X, T ) in Equation (1). The exponential decay in time associated to the fundamental solution is still evident in the case of fBm model: while it is difficult to notice it for the CTRW model. Studying the behavior of the fundamental solutions for V m for both the fBm and CTRW models in [Henry (2008)] it was determined that subdiffusive behavior enlarges the window of high potential at the soma, despite it lowers the maximum value of the peak, with respect to the classic solution in Equation (2), and it was suggested [Henry (2008)] that this effect could be in fact helpful to counterbalance deferred postsynaptic potentials over dendrites and to reduce temporal attenuation of the signal. Then, high spine density should have been related to more enhanced subdiffusive behavior.
More recently, other experiments have been performed on Purkinje cells and pyramidal cells [Santamaria(2011)] and the correlation between spine density and anomalous diffusion exponent in these types of neurons was explicitly studied; the anomalous MSD was described in terms of an exponent d ω by the introduction of a time dependent diffusion coefficient: and linear correlation was found between the parameter d ω and the measured density of spines in both the types of neurons studied. The effect of the system geometry on the transport regime in dendrites has been modeled considering the geometrical similarities between a comb structure and a spiny dendrite by the application of comb-like models of diffusion , ]. In this model, it was considered that particles may diffuse in both spines, the fingers of the comb, and the dendrite, the backbone of the comb, where spines behave as a traps for the moving particles, and the average survival time τ inside each spine is determined by its geometry. Markovian process was assumed inside each spine, i.e., exponential distribution of survival time Ψ M (t, τ ) = ∞ t 1 τ e −t ′ /τ dt ′ = e −t/τ , but the random size and shape of the spines [Nimchinsky(2002)] entail that the final process is the sum of many independent Markovian processes averaged over the distribution of the timescale τ , ]: When f S (τ ) is a power law Ψ(t) shows a power law behavior as well, and subdiffusive diffusion appears. In the present paper, we develop two models, CTRW and generalized grey Brownian motion (ggBm) like, to model diffusion in the system under study. Since the experiments described in the literature were performed using inert fluorescent tracers, the models proposed were kept as simple as possible and based on a diffusion process with leakage, to account for the loss of particles from the system. Subdiffusive behavior was included in both the models by means of the heterogeneity of the environment.

Results
The emergence of fractional kinetics in complex media in CTRW was introduced more explicitly as a general concept in [Pagnini(2014a)]. Analogously to the comb-like model presented, in that short note the special case of a survival probability of the Mittag-Leffer type was there derived in terms of a Markovian process with characteristic waiting time properly distributed: where E α (·) is the one parameter Mittag-Leffer function [Mainardi(2010)]: and f S (τ ) = 1 τ 2 K α (1/τ ). The distribution K α = K −α α,α is the fundamental solution of the space-time fractional diffusion equation [Mainardi(2001)] for the special case in which fractional orders of derivation of the space and time variables are equal, and of order α, with extremal asymmetry parameter equal to −α, leading to a solution defined on the positive real axes only. Within this approach f S (τ ) corresponds to the stationary distribution of these timescales. If instead we consider the non stationary case, it holds for the general case: however, in the non-stationary case the solution for f (τ, t) could be no unique given Ψ(t).
In the present case, the following identity holds: then for Ψ(t) = E α (−t α ) we may write: or equivalently: where M α (z) = W −α,1−α (z) is the M -Wright function, special case of the Wright function W λ,µ (z) defined by the series [Mainardi(2010)]: The relation in Equation (13) is a consequence of the Laplace transform relation between the M -Wright and the Mittag-Leffer functions [Mainardi(2010)]: thus: applying the change of variables q = rt α we have: Applying this idea to the most general solution for CTRW [Montroll(1964), Scalas(2004)] is it possible to write it in terms of a superposition of Markovian components, each characterized by the same jump PDF [DiTullio(2016)]: The simplest diffusion process of molecules associated to the transmembrane potential solution of the classic cable equation is P ′ M (r, t) = P M (r, q)e −q , with P M (r, q) = 1 √ 4πq e −r 2 /4q , standard diffusion process, multiplied by the exponential factor e −q that accounts for the loss of particles in the system. Following the same superposition principle after turning on the exponential decay term, the transmembrane potential P (r, t) corresponds to the integral of the solution of the classic cable equation averaged by the same H(q, t) = 1 t α M α (q/t α ), by considering that the decay is subjected to the same timescale of the diffusion process: The classic problem was written in terms of the adimensional variable T = t/τ , with τ = c m r m , and X = x/λ, with λ = r i /r m related to the circuit component of the membrane element. The solution of the fractional process can be written in terms of a superposition of the classic solution weighted by the distribution of the circuit element parameters, thus we have: Then in terms of circuit elements the system results characterized by a capacitance that varies between the elements of the circuit according to a certain time dependent distribution, considering r i , r m unitary constant for simplicity, in the present case it corresponds to: If there exists also a population of r m , representing the transmembrane resistance, the time decay of the solution and diffusion processes are described by two different but correlated distributions, because the coefficient r m disappears in the Gaussian factor.
In the comb-like model, the timescale is the average sojourn time in the spine and can be related to the geometry of the spine: a simple geometrical approximation involves spines composed by a head of volume V connected to the backbone by a neck of cylindrical shape with length L and radius a. The mean time spent in such a spine is τ = (LV )/(πa2D), where D represents a quantity called diffusivity of the spine , ]. If this volume may change dynamically it makes sense to consider a time-dependent distribution in this case as well.
Equation (21) can also be interpreted within the ggBm approach [Mura(2009)]; rewriting the integral form as follows: where inside the integral we recognize the fundamental solution expressed in Equation (6) of the fBm model defined in Equation (4) for the particular case α = κ. The ggBm-like stochastic process can be defined by the product: where X(t) is a Gaussian process with unitary coefficient of diffusion, rescaled by the diffusion coefficient D distributed according to: where comes natural the change of variables D = Λt α−1 , which is the fBm definition of the diffusion coefficient, thus p(Λ) = M α (Λ) [Molina (2016)]. The survival probability of each particle is conditioned to its diffusion coefficient D: The partial differential equation (PDE) for these processes can be derived by computing the Laplace-Fourier transform of the integral form in Equation (20), that readsV thus the transformed PDE is which correspond to the time fractional cable equation described in [Vitali(2017)] with 0 < α < 1 for Cauchy initial conditions: where ∂ α ∂T α is the Caputo time fractional derivative.

Discussion
Fractional calculus is often used to catch by parsimonious mathematical approach some underlying complex behavior. Caputo's fractional derivative is a non-local operator and for this reason, as pointed out in [Teka(2014)], it could be introduced to explain emergent behaviors such as the appearance of multiple timescale dynamics and memory effects, related to the complexity of the medium. In this work we derived two possible stochastic processes, CTRW and ggBm, for inert tracer diffusion in spiny dendrites that in principle give rise to the same partial differential equation for the transmembrane potential. The physical mechanisms expected behind CTRW and ggBm are strictly related to the process construction. Underlying the CTRW model there is the concept of trapping in the spines, where the distribution of timescales account the variability in the geometry and size of these spines. Underlying the ggBm model there is the idea that the environment is dynamical and that each particle may feel the surround in a different way. Both the approaches are approximations of the real system, and the aim of these models is to describe data behavior and possibly predict interesting biological features. The PDF evolution of both the processes is described by the time fractional generalization of the cable equation presented in Equation (29), that can be solved for the most common boundary and initial conditions by the application of the Efros theorem of Laplace transforms [Vitali(2017)]. The transmembrane potential function described by this model fulfills many of the biological features that have been previously suggested to explain spines role, as the compensation of delay in postsynaptic potentials and the attenuation time of the signal.
The first process is a CTRW built as a superposition of Markovian processes, each one subjected to a different timescale of the waiting time distribution, where the timescale follows a non stationary distribution. This means that the whole system change in time modifying the profile of the timescale distribution. The second process is based on a ggBm-like approach in which a Brownian process with unitary diffusion coefficient is rescaled by a random scale that is non stationary distributed as well. This scale represents in fact the value of the diffusion coefficient and can be used to generate anomalous time-scaling of the MSD in the final variable: If r(D, t) = 1 we obtain X ′2 ∼ t α . The exponential suppression has the effect that probability distribution collapse to zero in the infinite time, because all the particles disappear from the system.
Except for the exponential suppression that accounts for loss of particles, a similar non-stationary ggBm process has been proposed in [Molina(2016)] as an alternative to CTRW to account the Ergodicity Breaking (EB) described by several experiments on diffusion of cellular components in living systems. Despite both ggBm and CTRW may account for EB, in [Molina(2016)] it was shown that the p-variation test provides different values for the two alternative processes, and that values obtained for ggBm where compatible with the experimental dataset considered in their research, in contrast to CTRW.
For these reasons it seems promising to characterize the present processes looking forward for single particle tracking data to be compared with the models. Moreover, the two processes presented here account for the complexity of the phenomena directly from geometrical (waiting time timescales distribution) and/or electrophysiology (cell resistances and capacitance values distributions) properties of the system, that could be directly measured as it was done for spine density profiles in [Santamaria(2011)]. Finally, the anomalous transport phenomena is generated by a proper superposition of classic processes, that is not ad-hoc but can be related to experimental observations, clearly simplifying also the computational efforts of the simulation procedures of the trajectories.