Mathematical modeling of chemotaxis and glial scarring around implanted electrodes

It is well known that the implantation of electrodes for deep brain stimulation or microelectrode probes for the recording of neuronal activity is always accompanied by the response of the brain’s immune system leading to the formation of a glial scar around the implantation sites. The implantation of electrodes causes massive release of adenosine-5′-triphosphate (ATP) and different cytokines into the extracellular space and activates the microglia. The released ATP and the products of its hydrolysis, such as ADP and adenosine, become the main elements mediating chemotactic sensitivity and motility of microglial cells via subsequent activation of P2Y2,12 as well as A3A/A2A adenosine receptors. The size and density of an insulating sheath around the electrode, formed by microglial cells, are important criteria for the optimization of the signal-to-noise ratio during microelectrode recordings or parameters of electrical current delivered to the brain tissue. Here, we study a purinergic signaling pathway underlying the chemotactic motion of microglia towards implanted electrodes as well as the possible impact of an anti-inflammatory coating consisting of the interleukin-1 receptor antagonist. We present a model describing the formation of a stable aggregate around the electrode due to the joint chemo-attractive action of ATP and ADP and the mixed influence of extracellular adenosine. The bioactive coating is modeled as a source of chemo-repellent located near the electrode surface. The obtained analytical and numerical results allowed us to reveal the dependences of size and spatial location of the insulating sheath on the amount of released ATP and estimate the impact of immune suppressive coating on the scarring process.


Introduction
Nowadays, implantation of electrodes for deep brain recording of neuronal signals and electrical deep brain stimulation (DBS) are widely used for the therapy of several neurologic disorders [2][3][4]1]. In addition, increasing evidence suggests that implantable electrode technology, including multielectrode arrays (MEA), may be applied in the context of brain-machine interface [5][6][7]. However, recent studies reported a decrease in the quality of neuronal signals and the number of functional electrodes in chronic MEA recordings caused by the brain's immune response [8,9]. The in vivo response to the implantation of an electrode can be divided into the initial acute response to the cell damage and chronic response resulting in a formation of the electrode-tissue interface around the implantation site [8,10]. Several clinically used DBS devices use voltage-controlled stimulation, which leads to a strong dependence of the injected current on the impedance of the electrode-tissue interface, which separates electrodes from the surrounding brain tissue [14,16]. The impedance of the layer forming this interface increases with time, which in some cases makes it necessary to substantially increase the amplitude and/or frequency of the stimulation signal [18]. As follows from the immunohistochemical analysis, the main cellular components of the electrode-tissue interface are the microglial cells and reactive astrocytes [9]. Microglia are clustered around the implant in a reactive tissue sheath and persist throughout the whole lifetime on the amount of ATP released into the ECS. Apart from that, the impact of bioactive coating consisting of the IL-1 ra proteins on the glial scarring process is also taken into account in our model. The model of the glial scarring process is introduced in section 2. The analytical results for the linearized dimensionless equations are presented in section 3. In section 4, we analyzed the impact of the anti-inflammatory bioactive coating by means of stability theory. The results of numerical simulations for the two-dimensional case are presented in section 5. Our conclusions are given in section 6.

Model of microglial chemotactic acute response
The proposed macroscopic model of microglia chemotaxis is governed by the following equations:  3  4  3  2 2  2  2  1 1  1  5  3 3  3  3  2   AMP adenosine via cd73   4  5  4  3 3  3  3  2  4 4  4  4 adenosine uptake where n is the microglia density, y 1,2,3,4 are the concentrations of ATP, ADP, AMP and adenosine, respectively. D 1,2,3,4,5 are the diffusion coefficients and χ 1,2,3,4 are the chemotactic motility coefficients. We assume that the At the first stage, ATP released from the damaged cells (green filled circles) activates P2X 4,7 receptors and evokes an ATP-induced release of ATP from microglial cells [51] and thus results in an activation of P2Y 2,12 receptors, which provides signal amplification and facilitates the sensing of the shadow chemotactic gradients [27]. The second stage is characterized by the polarization of the cells within the chemotactic gradient field and autocrine activation of A3A receptors by adenosine formed from AMP by ecto-5′nucleotidase (CD73). At the last stage, adenosine (red filled circules) activates A2A receptors which causes retraction of the microglial cells and their movement towards the gradient of ATP.
chemotactic sensitivity of microglia to the ATP and ADP is concentration-dependent and is limited by the saturation of the binding process. The saturation of purinergic receptors is described by nonlinear interaction terms in the first equation and prevents the blow-up of a solution, which is observed in the classical Keller-Segel model [58]. The chain of kinetic equations describing the process of hydrolysis → y y ATP( ) ADP( ) 1 2 → → y y AMP( ) adenosine( ) 3 5 in the ECS [50] is padded by the equation for the concentration of IL-1 ra proteins y 5 . Here, the maximal rates of the enzymatic reactions are v 1,2,3,4 , w 1 , and k 1,2,3,4 , k i are the corresponding Michaelis-Menten constants. The maximal rate of ATP release by microglia is a. We use the simplifying assumption that amount of ATP released by the microglial cells is linearly proportional to the number of cells. Another additional source of ATP is the ATP-induced ATP release by reactive astrocytes whose maximal rate equals b [42]. The nonlinear function θ y ( ) 4 modulating sensitivity of microglia to adenosine is defined as the 2 . The specific form of the θ-function was chosen in such a way to take into account the interaction between adenosine and A3A/A2A receptors [27,34,54]. In particular, as can be seen from figure 2(a), the nonlinear θ-function takes positive values at moderate concentrations of adenosine in ECS, which means the enhancement of chemotaxis via activation of the adenosine A3A receptors, whereas rather high levels of adenosine cause an activation of A2A receptors which, in turn, leads to substantial retractions of microglial cells and their directed motion towards the electrode. The last equation of the system (1) describes the diffusion of IL-1 ra y 5 , which spreads away from the coated electrode surface and limits the amount of microglia aggregated around the implantation site. Apparently, the extent to which the IL-1 ra spreads in the ECS depends on the rate of its persistent emission L 0 and rate of interaction between microglia cells and IL-1 ra proteins d.

Analytical results for the linearized equations
Let us assume that the total number of microglial cells is fixed and consider the one-dimensional case with noflux boundary conditions:  where it is assumed that the electrode is placed at x = 0 and L is a finite interval or half-plane. We chose asymmetric boundary conditions for the ATP concentration: , where J 0 is the constant concentration of ATP, which is maintained by the persistent ATP release from microglia and astrocytes at the closest proximity to the electrode surface. Also, it is assumed that the concentration of IL-1 ra near the electrode surface is maintained at the constant level M 0 . Boundary conditions for ADP, AMP and adenosine are assumed to be zero . The initial conditions for the cell density and , and ξ t ( ) is the Gaussian noise of intensity ϵ.
The original model (1) can be non-dimensionalized by means of the following scales: and definitions: where = − n 10 3 is the equilibrium cell density. The time scale is normalized on a slow rate of adenosine uptake. The dimensionless equations are the following: ,ˆ,ˆ,ˆ,ˆ,     The studies of the relationship between the chemical signaling and nonlocal interactions of the cells clearly demonstrated that when the chemicals diffuse in the ECS more rapidly than the cells ( ≪ D 1 n ), the quasistationary spatial distributions of chemicals develop on a faster time scale [43,55]. The latter means stationarity of the spatial distributions of ATP, ADP, AMP, adenosine and IL-1 ra, which are described by the following steady states equations: Linear equations (7) can be solved one by one for the boundary conditions given above in (2). In the onedimensional case, the quasi-stationary spatial distributions of all chemicals y 1,2,3,4,5 can be written in the following explicit forms:    is the total gradient, maintaining the motion of microglia towards the source of the inflammation. As follows from the expressions (9) and (10), the density of microglial cells aggregated around the electrode is described by the following equation: The final expression for the cell density can be obtained after integration in (11): n As seen from figure 2(c), the instability of the homogeneous steady state leads to the appearance of an aggregate, the size and spatial location of which are determined by the coordinates of the global minimum and zero-crossing in grad(x) (10). The physical meaning of the zero-crossing in (10) can be easily understood as a point where the chemoattraction along the total gradient is balanced by the density of the glial scar. As can also be seen from figure 2(c), such zero-crossing points precisely correspond to the coordinates of the maximal values of the density in the formed aggregates. At the same time, the spatial coordinate of the global minimum in grad(x) might be considered as the point at which the influence of the cells' density overcomes chemoattraction and forms the natural 'backward border' of the aggregate. As follows from the results presented in figure 2(c), the usage of biologically active coat emitting IL-1 ra proteins leads to a drastic reduction of the microglial cells density inside the insulating sheath and shifts it away from the electrode surface. To study the impact of IL-1 ra on the microglia chemotaxis process, we calculated the size and absolute maximal values of the cell density inside the aggregate as functions of the amount of ATP, released from microglia and astrocytes.
As seen from figure 3(a), the spatial location of a global minimum x b defining the aggregate boundary strongly depends on the rate of the released ATP. At moderate rates of ATP release, the width of the insulating sheath is increased both in the presence of bioactive coating and when IL-1 ra is absent. The increase of the rate of ATP release in the absence of bioactive coating causes the growth and consequent saturation of the microglial cell density inside the aggregate (red solid line in figure 3(c)) as well as the shift of the point x max , where the cell density is maximal, in the direction towards the electrode (red solid line in figure 3(b)). The abrupt shrinkage of the sheath at some threshold value of J may be caused by drastic changes in grad(x) (10), meaning the growth of chemosensitivity with respect to ATP and ADP. The latter effect resembles, to some extent, the enhancement of the immune response, which is observed after the interaction of microglia and laminin [48] and is caused by the ATP-induced vesicle shedding and release of IL-1 cytokine by microglia [59]. In fact, the release of IL-1 cytokine by the microglia in response to the massive release of ATP might mean the appearance of another substance, playing the role of a chemo-attractant. It becomes clear in this context that the IL-1 ra electrode coating allows to weaken that part of the chemotactic signaling process, which is related to the ATP-induced release of IL-1 cytokine [59].
The application of the IL-1 ra coating causes qualitative changes in the tissue response to the growing level of extracellular ATP. In particular, the presence of bioactive coating weakens the additional chemo-attraction induced by the released IL-1 cytokines and prevents the shrinkage of the insulating sheath (green curve in figure 3(a)). Moreover, the presence of IL-1 ra induces strong reduction of microglial cell density inside the aggregate (green curve in figures 3(c)) and prevents its movement in the direction towards the electrode surface even for high levels of ATP release (green curve in figures 3(b)). Another important parameter affecting the process of microglia chemotaxis is γ, which is responsible for the functional regulation of adenosine A3A/A2A receptors. The substantial increase of the microglia cell density is observed for rather large values of γ in the absence of bioactive coating (see the red curve in figure 3(d)). This fact can be explained by the enhanced retractions of microglial cells which react this way to the increase of adenosine levels in the ECS. As shown in figure 3(d), the presence of IL-1 ra compensates the growth of the cell density caused by the increased level of extracellular adenosine. Dependences similar to those presented in figure 3 can be reproduced in a wide range of the values of parameters shown in table 1.
To study the role played by the parameters, characterizing the bioactive coating, on the aggregation process, we calculated the microglial cell density (13) for different values of the chemosensitivity coefficient χ 4 as well as for the different values of IL-1 ra a concentrations M 0 . As illustrated in figure 4, these two parameters have a similar effect on the cell density, namely, both a substantial increase of chemotactic sensitivity to IL-1 ra (see figure 4(a)) and a high concentration of IL-1 ra near the electrode surface lead to a strong reduction of the aggregate and limit the spread of inflammatory factors.

Stability of homogeneous steady state and impact of bioactive coating
To study the impact of the bioactive coating in terms of stability theory, we linearized the equations (6) about the homogeneous steady state n y y y y y (ˆ,ˆ,ˆ,ˆ,ˆ,ˆ) 1 (14) into (6) and linearization of the obtained system of equations allowed us to get the Jacobian matrix: 2 to obtain the sixth-order characteristic equation:  4(a)). Obviously, five eigenvalues take negative values for all wavenumbers. Thus, the spatial stability of the homogeneous steady state depends on the sign of λ 6 and is defined by the dispersion curve λ k ( ) 6 2 as shown in figure 5(a). To reveal the impact of the bioactive coating on the spatial stability of the homogeneous steady state, we calculated the dispersion curves λ k ( ) 6 2 for different values of the chemosensitivity χ 4 . As follows from the results presented in figure 5(b), an increase of the chemosensitivity to the IL-1 ra leads to a substantial decrease of the area of instability where λ 6 takes positive values. Moreover, the analysis of the roots of the characteristic equation (15) by means of the Routh-Hurwitz theory allowed us to estimate the boundaries of the domain of the spatial instability. The results of this analysis are in good agreement with the results of the direct numerical solution of equation (15). As illustrated in figure 5(b), the area of instability is finite, localized in space near the electrode and decreases as the coefficient of chemosensitivity to the IL-1 ra χ 4 increases. Also the absolute values of λ 6 decrease in the presence of the IL-1 antagonist, what may be interpreted as a reduction of the instability caused by the  release of ATP in ECS. These last findings are in good agreement with the analytical results presented above concerning the impact of the bioactive coating on the size and location of the glial sheath. Indeed, as seen from the results presented above in figure 2(c), the use of the bioactive coating leads to a substantial decrease of the cell density inside the aggregate as well as its removal from the electrode.

Results of numerical simulations
To verify the analytical results obtained for the original nonlinear model, we performed numerical integrations of the dimensionless equations (5) for the two-dimensional case by using the standard finite-difference scheme together with the high order explicit method for stiff differential equations utilizing the explicit Runge-Kutta-Chebyshev formulas [52]. The choice of the two-dimensional geometry may simplify the comparison of the obtained numerical results and results of representative microscopy (X10) data of horizontal brain sections immunostained for GFAP surrounding the electrode insertion site, which are presented in [49]. We assumed that the electrode is located in the center of the square domain which allows us to model the realistic situation, when it is symmetrically enwrapped by the microglial cells. The obtained values of the cell density were normalized by the total density inside the chosen domain, which is calculated as a sum of all values of n, calculated in each point on the equidistant grid. Since the expressions for A 1 , A 2 , A 3 are different for the systems (5) and (6), then for the numerical integration of the system (5)  , respectively. As shown in figure 5, the emission of IL-1 ra proteins by the anti-inflammatory coating leads to a drastic reduction of the cell density around the electrode and nearly prevents the scarring process. The two-dimensional stationary spatial distributions of ATP, ADP and adenosine have similar profiles to those obtained in the onedimensional case shown in figure 2(b). These findings are in good qualitative agreement with recent experimental results [49]. In particular, the two-dimensional plots in figure 6 resemble the results of representative microscopy (X10) data of horizontal brain sections, which can be found in [49]. As already discussed in section 3, the strategy of the suppression of the acute immune response aims to eliminate the chemotactic signaling which is related to the release of IL-1 cytokines and may be considered as a very effective approach to prevent scar formation. Moreover, the suppression of the cytokine release may play a crucial role in the prevention of astrogliosis responsible for substantial decrease of the extracellular volume. The shrinkage of extracellular volume is known to induce the appearance of diffusive barriers in the ECS which are able to affect the diffusion of molecules and limit the spread of the chemotactic signaling to the closest proximity of the electrode [60].

Conclusions
We presented a simple mathematical model of the acute chemotactic response of the microglia to the electrode implantation, which is accompanied by a massive release of ATP. Given the recent experimental results, we assumed that the chemotactic signaling and motility of microglia are maintained by the chain of enzymes' reactions on the surfaces of microglial cells, which performs the hydrolysis of ATP to adenosine. Based on our model, we found that the release of a substantial amount of ATP into the ECS induces the chemotactic motion of microglia towards the electrode. The stationary spatial distributions of ATP and its hydrolysis products together with the stationary distribution of microglial cell density were calculated analytically for the one-dimensional case. We also demonstrated the possibility of an abrupt shrinkage of the insulating sheath for some threshold values of ATP concentration.
With our model, we analyzed possible effects of anti-inflammatory bioactive coating on the glial scarring process. In particular, we found that a recently proposed bioactive coating shifts the emerging encapsulating sheath away from the electrode, preventing its abrupt shrinkage and stabilizing the cell density around the electrode at a moderate level. Thus, the usage of the anti-inflammatory bioactive coating containing the proteins of the IL-1 antagonist, is able to suppress or significantly diminish the scarring process around the implantation site which, in turn, can increase the conductance and signal-to-noise ratio of the electrode-tissue interface.