{'fraction_non_alphanumeric': 0.02831705690865592, 'fraction_numerical': 0.008706906694905283, 'mean_word_length': 4.902314814814815, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A model of fractional-order diffusion in the glial scar

Implantation of neuroprosthetic electrodes induces a stereotypical state of neuroinflammation, which is thought to be detrimental for the neurons surrounding the electrode. Mechanisms of this type of neuroinflammation are still not understood well. Recent experimental and theoretical results point out possible role of the diffusion species in this process. The paper considers a model of anomalous diffusion occurring in the glial scar around a chronic implant in two simple geometries -- a separable rectilinear electrode and a cylindrical electrode, which are solvable exactly. We describe a hypothetical extended source of diffusing species and study its concentration profile in steady-state conditions. Diffusion transport is assumed to obey a fractional-order Fick law, which is derived from physically realistic assumptions using a fractional calculus approach. The derived fractional-order distribution morphs into regular order diffusion in the case of integer fractional exponents. The model presented here demonstrates that accumulation of diffusing species can occur and the scar properties (i.e. tortuosity, fractional order, scar thickness) can influence such accumulation. The observed shape of the concentration profile corresponds qualitatively with GFAP profiles reported in the literature. The main difference with respect to the previous studies is the explicit incorporation of the apparatus of fractional calculus without assumption of an ad hoc tortuosity parameter. Intended application of the approach is the study of diffusing substances in the glial scar after implantation of neural prostheses, although the approach can be adapted to other studies of diffusion in biological tissues, for example of biomolecules or small drug molecules.


Introduction
Implantation of neuroprosthetic electrodes induces a sustained state of neuroinflammation and scarring, which is thought to be detrimental for the neurons surrounding the electrode [18]. Literature demonstrates that in chronic conditions the recording longevity of such electrodes in experimental animals is highly variable (i.e. for wire electrodes - [13], silicon-based electrodes and multiwire arrays - [38]). Over 100 studies have described stereotypic features of the brain response to microelectrodes that occur irrespective of the type of implant, method of sterilization, species studied, or implantation method [8].
The formation of the glial scar is a complex reactive process involving interactions between several types of cells, notably astrocytes and activated microglia, which are mediated by plethora of bio-active molecules (i.e. cytokins). The reactive astrocytes form a dense web of interdigitated processes which over-expresses the Glial Fibrillary Acidic Protein (GFAP) and attach to the implant (see Fig.  1). Thus GFAP is commonly used in neuroprosthetic studies as a marker of neuroniflammation. In parallel, the microglia are rapidly activated in a wide area around the lesion site and undergoe a profound change in cell shape and phenotype. Concurrent with the glial scar formation, neuronal density within the recording radius of the microelectrodes is reduced, leading to even fewer distinguishable single-unit recordings [5,37,10,1,31]. During this process the cells change substantially the composition, the morphology and the functional properties of the extracellular matrix. Roitbak and Syková [33] also demonstrate changes of the diffusion path in reactive astrogliosis states and therefore in the extracellular space (ECS) properties.

Figure 1: GFAP distribution around an implant site
A -GFAP staining after 6 weeks of implantation in tethered configuration; B -Mean intensity distribution as a function of the distance to the insertion track [30]. The dataset was published previously in [39]. Scale bar -200 um.
Several hypotheses for these observations have been put forward. Some authors proposed that formation of glial scar results in a physical barrier to diffusing substances, thus creating a toxic environment for the neurons. Others have proposed that many neurons around the electrodes die shortly after implantation [5,1]. More recently, McConnell et al. [18] have proposed that the observed loss of signal can also result from the progressive degeneration of nerve fibers and synapses due to persistent local chronic inflammation. Ward et al. attribute electrode failure to the traumatic injury resulting from insertion and a long-term foreign body response to the implant [38].
While there are observations supporting each hypothesis there is no agreement as to which is the predominant effect at different time scales (review in [26,12]). Quantitative histological descriptions of the glial scar may therefore demonstrate these effects more convincingly and provide metrics for more robust safety and biocompatibility assays for neural prostheses.
In this paper we study accumulation of a diffusing species in steady-state conditions, produced as a reaction to the presence of another object, i.e. an implanted electrode. We hypothesize that such conceptual model can describe diffusion phenomena occurring in the extracellular matrix following implantation of flat electrodes, for example Michigan type of silicon probes or cylindrical electrodes, such as microwires.
2 Extended source compartment model Implanted electrodes have simple cross-section geometries and high aspect ratio with regard to height. In such way, implant geometry can be approximated as consisting of ideal shapes, i.e. an infinite cylinder, a line or a half-plane. This allows for applying a symbolic calculation approach in the modeling of the diffusion problem.
As an idealized situation we will consider semi-infinite medium with two spatial compartments: a source compartment (S) and a tissue compartment (T), (see Fig. 2). To simplify calculation we will assume that the S-compartment will have only constant production of diffusing species but no degradation, while the T-compartment will have only first order degradation but no production. The overall solution will be represented by a sum of two terms for the two coupled domains. In the following presentation the term c s will represent the concentration of the substance in in the glial scar, i.e. the S-compartment. while the second term c f will represent the concentration in the outer zone -i.e. the T-compartment.
Further, the terms are assumed to be orthogonal, i.e. c s (x) = 0 in the distal compartment x > L whereas c f (x) = 0 in the proximal compartment x < L.
Since the overall solution is assumed to be continuous the following conditions will have to be imposed on the boundary between domains: In addition, we also assume natural boundary conditions in steady statec(0) = c 0 , c(∞) = 0. Subsequent analysis is performed both in traditional partial differential equation setting and in fractional calculus (i.e. differ-integral) approach as a possible conceptual generalization.

Anomalous diffusion in complex media
At present it is well established that the ECS occupies a volume fraction of between 15 and 30% in normal adult brain tissue with a typical value of 20% and that this falls to 5% during global ischemia [36]. Therefore, the diffusion impediment can not be neglected in modeling of biological diffusion. Impeded diffusion can be modeled by two classes of models.

Regular impeded diffusion
In this framework, the medium imposes only a spatial impediment on the diffusing species and does not change the diffusion law. Nicholson et al [22,21] consider that the densely packed cells of the brain and their interstitial spaces resemble a porous medium with two phases, an intra-and extracellular phase. Diffusion in the permeable phase of porous media is analogous to the diffusion in the narrow spaces between brain cells, i.e. the extracellular space, ECS. In accordance, volume transport of species having concentration c is governed by the rescaled reaction-diffusion equation: where the source term is denoted by s (1/concentration) , the clearance term denoted by k(c) = κ c is assumed to be first order decay and D is the diffusion coefficient (length 2 / time). The following additional parameters are represented respectively by the parameters λ -tortuosity (dimensionless); α -porosity coefficient (dimensionless); κ -clearance speed ( 1/time).
Unless stated explicitly, we will further use a re-parametrization of the problem preserving the form and interpretation of the equations:

Fractional diffusion phenomena
Since the brain extracellular space is a complex medium this transport equation can be regarded only as an linear approximation. In accordance, tortuosity can be considered as the linear correction for the anomalous diffusion [21] and [36]. Diffusion in porous media, such as tissues, are characterized by deviation from the usual Fick's diffusion laws. Notably, these processes are characterized by a distribution of waiting times, having heavy tails approximated by an inverse fractional power law, and hence, dependence of the effective diffusion phenomenon on the duration of the measurement. Also the law of the mean squared displacements exhibits deviations from linearity. The resulting behavior can not be described by an ordinary differential equation.
Since the end of the XX th century a different approach to handle deviations from the classical diffusion and hence the anomalous diffusion phenomena has been employed. In this framework the underlying physical process is modeled as a continuous time random walks [20], using the mathematical apparatus of fractional calculus [24]. In the fractional calculus approach the transport equation reads where ∇ β denotes the fractional-order Fick's law, i.e. the fractional flux in the system. Interpretation for the symbol of the last operator will be given further in the relevant sections of the paper. This form of the fractional Fick's law naturally implies spatial non-locality and can be derived from rigorous approaches using spatial averaging theorems [19,40].

Regular diffusion along the line
In this case we consider the following system: In steady state the equations have the following general solutions: Then since the solution is limited k 1 = 0. After some algebraic transformations we arrive at the following algebraic system to be solved The maximum concentration is attained at It is easy to check that at x = 2 x m we have c s (2 x m ) = c 0 .

Fractional diffusion along the line
The fractional Fick's law in this case reduces to the fractional derivative of a function in the sense of Caputo [3], [4]. This derivative is in turn defined by the differ-integral: where Γ(.) is the Euler's gamma function, which for integer number is equal to the factorial Γ(n) = (n − 1)! The fractional Fick's law is given in this case by For the convenience of the reader a simple derivation is presented in Appendix A. Combining this equation with a conservation of mass equation for the concentration of particles leads to the fractional diffusion equation of the type.
For the T compartment since the first derivative of the equation is constrained at the boundary we will reformulate the problem entirely in terms of Caputo derivatives.
The equation can be solved in terms of special functions [17,16]. The general solution is given by where the E a,b denotes the Mittag-Leffler function (see Appendix B ). To avoid unbounded solutions at infinity we pick up constants with opposite signs. Therefore, finally where the constant C can be determined from the boundary or initial conditions. Since κ D is assumed to be small we will use further the asymptotic expression valid for small values of x to illustrate the relationship to the integer-order solution. c The resulting system can be solved by subsequent integration using eq. 6 to yield the following algebraic system for the first boundary and second boundary condition, respectively. In the resulting system k 2 and k 5 are unknown constants to be determined by the initial and boundary conditions. Considering further the steady state gives finally for the proximal compartment while for the distal component 5 Cylindrical geometry

Regular cylindrical diffusion
In the case of implanted wire electrodes, we will assume cylindrical symmetry of the problem. The the Laplacian operator can be represented as We will look only for solutions that are cylindrically and axially symmetric that is ∂c ∂φ = 0 and ∂c ∂z = 0. Therefore, we have to consider the following system: The solution will be presented as a sum of two terms for the two coupled domains where we also impose continuity conditions at the border: In addition we assume also a natural boundary conditions in steady state: In steady state the equations transform to the following independent system The general solution is given by From the finiteness of the solution at the left boundary it follows that c 2 = 0 and c s (r) = c 0 − s 4D r 2 (1)

Compartment T The equation can be transformed in the form
which corresponds to the modified Bessel equation [27] by the ansatz q 2 c f (qr) = y(r). Notably, The general solution of this equation is given in terms of the modified Bessel functions of the first and second kind, notably : Since I 0 (r) diverges at infinity the only acceptable solution is Applying the continuity conditions results in the system The system can be solved to yield: From the explicit form of solution it can be seen that the value of the concentration at the origin is completely fixed by the geometry of the problem.

Fractional cylindrical diffusion
Generalization to the 2D case is a considerably more difficult problem. Following the derivation of Meerschaert [19] we define the fractional Fick's law of fractional order β < 1 as where J β denotes the vector fractional integral operator.
Compartment S Following the derivation of Meerschaert [19] we define the fractional Fick's law. Therefore, in steady state we have D ∇ · ∇ β c s + s = 0 The equation can be solved by Laplace's transform method to yield Similar considerations about the boundedness of the solution at r = 0 require that c 1 = 0.
Compartment T In steady state conditions we have The Laplace transform method can be also applied in this case, although no explicit solution can be identified ( see Appendix C.2). Nevertheless, for completeness of the presentation a similar procedure is applied. This is equivalent to the assumption that the fractional exponent equals two in the T domain.
Applying the continuity conditions results in Let's assume that the electrode is thick, with a finite radius ρ. This case is known as the standard Neumann boundary condition ∂ c s (r) ∂r r=ρ = 0 Linear geometry In the case of tight contact between the tissue and the electrode one can assume as well vanishing of the flux at the boundary of the electrode. In this case we will have which gives immediately for the concentration at the boundary The maximum of the last expression is attained at x = ρ.

Cylindrical geometry The general solution for the S-compartment is given by
After some algebraic manipulations the solution can be simplified to It can be shown that in the limit ρ → 0 we can recover In the case of fractional diffusion after some algebra we can obtain which again for ρ → 0 and β → 1 reduces to Eq. 1.

Non-vanishing flow
Let's denote by z the value of the flow attained at the electrode boundary. In this case we have to consider the condition ∂ c s (r) ∂r r=ρ = z Linear geometry In this case we have the condition After some algebra we arrive at: The position of the peak is attained by solving the equation which gives the unique value Cylindrical geometry In this case we have the condition and after some algebraic manipulations It can be shown that in limit ρ → 0 we can recover The position of the peak is attained by solving the equation which gives the physically meaningful value of From this analysis it is apparent that the peak position is influenced by the direction of the flux.

Numerical analysis
Some numerical results are demonstrated in Fig. 3 using literature values of parameters for rat neocortex, [21]. Available experimental data show that a typical value for small molecules for λ in the nervous system is about 1.6, which implies thatD is some 2.6 times smaller than D. For macromolecules λ is increased, in part because of more frequent interaction with the walls of the narrow channels through the ECS [23]. Numerical values used in the calculation are as follows: D = 1.10 −6 cm 2 s −1 for proteins, λ = 2.5 (NGF) and α = 0.21; uptake k = 1.10 −4 (NGF). Magin et al. [15] measure experimentally the fractional order of 1.95 using diffusion MRI. This corresponds to the value β = 0.95 in our parametrization.
All cases use the same values of porosity. The glial scar thickness was taken as the range L = 50−150 µm. For the unconstrained diffusion case no correction for λ is given. The linearized case corresponds to the approach of Nicholson and Sykova [21,36]. In the fractional cases no tortuosity correction is applied but only the fractional order β is varied. Results are plotted in Figs. 3 and 4.

Diffusion around brain implants
Presented results show that in steady-state conditions substances, which are produced presumably in the glial scar could accumulate and produce characteristic profiles, matching qualitatively empirically observed GFAP and microglial profiles published in available experimental literature. More specifically, shapes of the concentration profiles correspond qualitatively with the GFAP profiles, for example [44,18,42,41,29,28,35], and the blood vessel area distribution [6] reported in literature. In all of these studies implanted probes had planar geometry, therefore reported GFAP profiles can be discussed under the assumptions of the presented model. Presented model can be extrapolated to separable 3D geometries as well considering presence of planar symmetries. The shape of the cylindrical distribution corresponds also to studies employing cylindric or wire electrodes [32,10]. Presented approach, therefore, provides means for approximation of such profiles and estimation of geometrical parameters related to the problem, such as the scar thickness and the intensity of the source. In some cases the fractional exponent can supplement the so-far used linearized parameters, such as tortuosity, if such estimates are not available.
Our model demonstrates that ex-centric accumulation of diffusing species can occur and the scar geometry (i.e. tortuosity, fractional order, scar thickness) can influence such accumulation. From the analysis of the boundary condition we can conclude that the insulation of the electrode with regard to the boundary (i.e. presence or absence of lateral flux) determines the shape of the observed profile in both studied geometries. We can therefore speculate, that all studies exhibiting profiles with ex-centric accumulation actually had somehow loose contact between the implant and the brain tissue, possibly due to micromotion. To further verify that, however, additional efforts in mechanical modeling of devices are necessary.
Recently [35] modeled numerically diffusion around implants in order to design a diffusion sink placed at the device surface that would retain proinflammatory cytokines for sufficient time to passively antagonize their impact on the foreign body response. Such setting would correspond to a negative flux condition in our model and it can indeed be shown that this will result in reduction of the gradient. Authors used no-flux boundary conditions and presented results, which qualitatively correspond to the analysis presented here. On the other hand, reported experimental results demonstrate deviations from the numerical model, which can be interpreted well using our model.

Impeded diffusion phenomena
Anomalous diffusion problems naturally arise in the setting of complex biological environment. Modeling of diffusion in different complex media could provide further understanding in a variety of experimental conditions. Impeded or anomalous diffusion in the brain is already a well established phenomenon [36]. Sykova and Nicholson [36] list several factors that determine diffusion impediments: i) an increase in geometric path length; ii) transient trapping of molecules in dead-space microdomains; iii) an increased interstitial viscous drag on migrating molecules; iv) transient binding to membrane-attached or extracellular matrix-attached receptors; and v) nonspecific interaction with fixed opposite charges. All five factors can be thought of as introducing a delay into the passage of a molecule in brain tissue relative to that in a free medium. Therefore, this complex medium can be conceptually modeled as the presence of two phases -one permeable and one impermeable, which impedes undergoing diffusion.
Fractional diffusion models have been employed in hydrology describing well slow diffusion [19] in protein diffusion in the plasma membrane [11,9]. The spatial complexity of a medium can impose geometrical constraints on transport processes on all length scales that can fundamentally alter the laws of standard diffusion [20]. Specifically fractional order diffusion models have been employed to model water diffusion in the brain [14].
Fractional models can provide new insights on mesoscopic aspects of the studied phenomena. For example, space fractional models have been used to describe cardiac cell conduction [2], and diffusion phenomena in fixed tissue samples [15].
ECS of the brain comprises the matrix that resides outside the neurons and glia cells, most importantly the interstitial space between neighboring cells. The ECS is a reservoir for ions involved in electrical activity, a communication channel for chemical messengers and a conduit for drug delivery. A quantitative description of extracellular diffusion is important whenever the transport of neurotransmitters, neuromodulators, and therapeutics in the brain ECS is considered. Diffusion-mediated transport of biomolecules substances is hindered by the ECS structure but the microscopic basis of this hindrance is not fully understood [7]. Evidence for anomalous diffusion in the brain has been found in the rat cerebellum [43], which provides support also for employing fractional models in the continuous approximation limit. Most probably, dead space microdomains can be the cause of such anomalous diffusion [34].
The magnitude of the exponent β can be used as a metrics for the departure from linearity. Since decrease of β will result in behavior approximating advection flow, this can used to reason about state where such flows can be enhance, for example in the case of the leaky blood-brain barrier or in conditions with increased local blood flow.
Both definitions coincide for problems where the function and its first n derivatives vanish at the lower limit of integration, i.e. when f (a) = 0, . . . , f (n) (a) = 0. Caputo's definition is better suited for problems where the function and its derivatives doe not vanish at this limit, because in this case it is given as kind of regularization of the Riemann-Liouville differintegral to avoid divergences.
Caputo's derivative is a left inverse of the fractional integral: while the fractional integral is a conditional inverse of Caputo's derivative: in the case when D β a f = 0. This allows one to solve simple fractional systems, such as the fractional diffusion equation in the example.
Fractional Fick's law for an exponent β > 0 can be derived under the assumption that the transport is given by a wighted fractional average of the differential of concentrations in the domain [a, x]. To avoid unphysical divergence at the borders of the domain, which are of interest in our problem we further regularize by subtracting the boundary concentration c(a): One can recognize this expression as the usual definition of the Riemann-Liouville fractional derivative. The integral can be evaluated partially assuming as usual existence of the derivative of the concentration c (t). Then applying integration by parts we get where we notice that the first term evaluates to zero.
Differentiating the last integral by x gives However, the last expression can be recognized as D 1−β a c. Therefore,

B Mittag-Leffler functions
During the last decades fractional calculus has been widely applied in many scientific areas ranging from mathematics and physics, up to biology, engineering, and earth sciences. The Mittag-Leffler functions play an important role in fractional calculus since many solutions of fractional differ-integrals can be expressed in terms of Mittag-Leffler functions. The one parameter ML function is a generalization of the exponential function: The two parameter ML function is an additional generalization in the sense , α > 0 A particular form of the one parameter ML function studied by F. Mainardi is e α (t) = E α (−t α ) The following result can be stated D α 0 e α (b t) = b e α (b t). It is common to point out that the function e(t) matches for t ≈ 0 a stretched whereas as t → ∞ with a negative power law. The short time approximation is derived from the power series representation as follows [16] e α (t) ∼ exp − t α Γ(1 + α) , t → 0 e α (t) ∼ Transforming back to the spatial domain this gives a general solution in terms of modified Bessel functions: f (r) = c 2 I 0 (r) + c 1 K 0 (r) .