Asymptotic and numerical analysis of a stochastic PDE model of volume transmission

TL;DR

This paper develops an asymptotic and numerical analysis of a stochastic PDE model for volume transmission, revealing that neurotransmitter concentration is primarily spatially uniform and independent of many details, with second-order effects influenced by various factors.

Contribution

The paper provides the first analytical asymptotic expansion for a stochastic PDE model of volume transmission, highlighting the spatial uniformity of neurotransmitter concentration and its dependence on second-order factors.

Findings

Concentration is spatially constant to leading order.

Leading order concentration is independent of nerve varicosities and geometry.

Second order effects depend on firing correlations and space geometry.

Abstract

Volume transmission is an important neural communication pathway in which neurons in one brain region influence the neurotransmitter concentration in the extracellular space of a distant brain region. In this paper, we apply asymptotic analysis to a stochastic partial differential equation model of volume transmission to calculate the neurotransmitter concentration in the extracellular space. Our model involves the diffusion equation in a three-dimensional domain with interior holes that randomly switch between being either sources or sinks. These holes model nerve varicosities that alternate between releasing and absorbing neurotransmitter, according to when they fire action potentials. In the case that the holes are small, we compute analytically the first two nonzero terms in an asymptotic expansion of the average neurotransmitter concentration. The first term shows that the concentration is spatially constant to leading order and that this constant is independent of many details in the problem. Specifically, this constant first term is independent of the number and location of nerve varicosities, neural firing correlations, and the size and geometry of the extracellular space. The second term shows how these factors affect the concentration at second order. Interestingly, the second term is also spatially constant under some mild assumptions. We verify our asymptotic results by high-order numerical simulation using radial basis function-generated finite differences.

Equations

$$\frac{\partial c}{\partial t}$$

$$c(x,t)$$

$$D{\partial }_{\nu }c(x,t)$$

$${\partial }_{\nu }c(x,t)=0,\text{if }x\in \partial E\backslash {\cup }_{n=1}^N\partial V_n.$$

$$\underset{t\to \infty }{\lim }\mathbb{E}[c(x,t)]$$

$${\theta }^{(1)}$$

$${\theta }^{(1)}=\frac{{\mu }_D}{{\mu }_C},$$

$${\mu }_D=\frac{1}{N}\sum _{n=1}^N{D}_n{\pi }_1^{(n)},{\mu }_C=\frac{1}{N}\sum _{n=1}^N{C}_n{\pi }_0^{(n)}.$$

$$|x_n-x_m|=\mathcal{O}(1)\text{for }n\ne m,$$

$$|x_n-x|=\mathcal{O}(1)\text{for }n\in \{1,\dots ,N\}\text{ and }x\in \partial U.$$

$$E:=U\backslash \overline{{\cup }_{n=1}^N{V}_n}.$$

$$J(t)=(J_1(t),\dots ,J_N(t))\in \mathcal{J}\subseteq \{0,1{\}}^N,$$

$$\mathcal{J}=\{(0,\dots ,0),(1,\dots ,1)\}⊊\{0,1{\}}^N.$$

$${\rho }_{\mathbf{i}}{Q}_{\mathbf{i},\mathbf{j}}={\rho }_{\mathbf{j}}{Q}_{\mathbf{j},\mathbf{i}},\text{for all }\mathbf{i},\mathbf{j}\in \mathcal{J}.$$

$$⟨\mu ,\eta ⟩:=\sum _{\mathbf{i}\in \mathcal{J}}{\rho }_{\mathbf{i}}{\mu }_{\mathbf{i}}{\eta }_{\mathbf{i}},\mu ,\eta \in {\mathbb{R}}^{|\mathcal{J}|}.$$

$$\frac{\partial c}{\partial t}$$

$${\partial }_{\nu }c$$

$$c$$

$${\partial }_{\nu }c$$

$$v_{\mathbf{j}}(x):=\underset{t\to \infty }{\lim }\mathbb{E}[c(x,t)1_{J(t)=\mathbf{j}}].$$

$$1_{J(t)=\mathbf{j}}=\{\begin{array}{ll}1& \text{if }J(t)=\mathbf{j},\\ 0& \text{if }J(t)\ne \mathbf{j}.\end{array}$$

$$\Delta \mathbf{v}(x)+Q^T\mathbf{v}(x)=0,x\in E,$$

$${\partial }_{\nu }\mathbf{v}$$

$$\begin{array}{rl}{v}_{\mathbf{j}}& =0,x\in \partial V_n\text{if }{j}_n=0,\\ {\partial }_{\nu }{v}_{\mathbf{j}}& ={\rho }_{\mathbf{j}},x\in \partial V_n\text{if }{j}_n=1,\end{array}$$

$$(I-B_n)\mathbf{v}+B_n({\partial }_{\nu }\mathbf{v}-\rho )$$

$$(B_n{)}_{\mathbf{i},\mathbf{j}}=0\text{ if }\mathbf{i}\ne \mathbf{j},\text{and}(B_n{)}_{\mathbf{j},\mathbf{j}}=j_n\in \{0,1\},\text{where }\mathbf{j}=(j_1,\dots ,j_N).$$

$$P^{-1}{Q}^{T}P=-\Lambda =-\text{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{|\mathcal{J}|}),$$

$$0={\lambda }_1<{\lambda }_2\le \cdots \le {\lambda }_{|\mathcal{J}|},$$

$$\mathbf{u}:=P^{-1}\mathbf{v},$$

$$\Delta \mathbf{u}-\Lambda \mathbf{u}=0,x\in E,$$

$${\partial }_{\nu }\mathbf{u}=0,x\in \partial U,$$

$$(I-B_n)P\mathbf{u}+B_n(P{\partial }_{\nu }\mathbf{u}-\rho )$$

Full text

{AMS}

35R60, 60H15, 35B25, 35C20, 65M60, 65D25,

1 Introduction

Synaptic transmission (ST) is a fundamental neural communication mechanism. In ST, a neuron fires an action potential that travels down its axon to a synapse that is adjacent to a second neuron. Upon arrival at the synapse, this electrical signal is translated into a biochemical response that results in neurotransmitters being released into the synaptic cleft. The neurotransmitters then rapidly diffuse across the synaptic cleft and bind to receptors on the second neuron, which then makes the second neuron more or less likely to fire, depending on whether the neurotransmitters are excitatory or inhibitory. Notice that ST is a “one-to-one” mode of communication, as it occurs through an essentially private channel between two neighboring neurons [1].

Complementing ST, there is an important “many-to-many” neural communication pathway called volume transmission (VT) [2, 3, 4]. In VT, a collection of neurons projects to a distant brain volume (a nucleus or part of a nucleus) and increases the neurotransmitter concentration in the extracellular space in the distant volume by firing action potentials [4]. The effect of this process is the modulation of ST in the distant volume, and for this reason VT is sometimes referred to as neuromodulation. Two important examples of VT are the serotonin projection from the dorsal raphe nucleus to the striatum [5, 6] and the dopamine projection from the substantia nigra to the striatum [7] (see [4] for more examples). VT is thought to play an important role in maintaining the sleep-wakefulness cycle, motor control, and treating Parkinson’s disease and depression [4]. Furthermore, very recent experiments suggest that VT may be a common neurobiological mechanism for the control of feeding and other behaviors [8].

In spite of the biological and medical significance of VT, mathematical modeling has only recently been used to study this process [9, 10, 11]. In contrast, mathematical modeling has been used to study ST for decades. Indeed, the model of Hodgkin and Huxley transformed the field of neuroscience (Hodgkin and Huxley won the 1963 Nobel Prize in Physiology or Medicine) and continues to influence current research [12, 13]. Furthermore, modeling ST has led to innovations in dynamical systems theory [14] and still prompts new mathematical questions [15, 16, 17].

In this paper, we use a stochastic partial differential equation (PDE) model of VT to investigate the average neurotransmitter concentration in the extracellular space. In particular, we are interested in how the average neurotransmitter concentration varies across the space and how it depends on (i) the number, arrangement, and shape of nerve varicosities, (ii) neural firing statistics, (iii) the amount of neurotransmitter released when a neuron fires, (iv) the neurotransmitter diffusivity, and (v) the size and geometry of the extracellular space. We answer these questions in the case that the nerve varicosities are small compared to the distance between them. We discuss the biological relevance of our model in the Discussion section.

Mathematically, we use the diffusion equation in a bounded three-dimensional domain to model the neurotransmitter concentration in the extracellular space. This domain has a set of interior holes that represent nerve varicosities, see Figure 1. Varicosities are bulbous enlargements on an axon that contain neurotransmitter [18]. When a neuron fires, neurotransmitter is released from a nerve varicosity into the extracellular space, and when a neuron is not firing (quiescent), neurotransmitter is taken back up into the varicosity from the extracellular space. Hence, each of these interior holes has boundary conditions (either an inhomogeneous flux condition or an absorbing condition) that switch according to neural firing, which we assume follows a continuous-time Markov jump process.

We first apply matched asymptotic analysis to a large system of deterministic PDEs that describes the steady-state mean of the stochastic PDE. Our asymptotic analysis yields the first two nonzero terms in an asymptotic expansion of the average neurotransmitter concentration, where the small parameter compares the size of varicosities to the distance between them. The first term in this expansion is constant in space, and this constant depends on very few details in the problem. Interestingly, if all the neurons have the same size, shape, and proportion of time firing, then the second term in the expansion is also constant in space. We verify our asymptotics by comparing to a simple example which can be solved analytically and by comparing to detailed numerical solutions for more complicated examples.

We now comment on how the present study relates to prior work. A similar stochastic PDE model of VT in one space dimension was proposed and analyzed in [11]. This model was generalized to two and three space dimensions in [10]. In [10], probabilistic methods were used to prove that the average neurotransmitter concentration is constant in space to leading order for small varicosities and to compute this constant. The asymptotic methods of the present study allow us to compute the next order term for the average neurotransmitter concentration and to relax the assumption in [10] that the varicosities are spherical. In further comparison to [10], in the present study we develop numerical methods to study this stochastic PDE model.

More broadly, the present study employs and develops techniques from three disparate areas of mathematics. First, we use a recently developed theory that yields deterministic boundary value problems that describe statistics of parabolic PDEs with stochastically switching boundary conditions [19, 20, 21]. This theory has been used in diverse areas of biology, including insect respiration [22], intracellular reaction kinetics [23], and transport through gap junctions [24, 25, 26]. Second, our asymptotic analysis uses strong localized perturbation theory [27, 28, 29], which applies to perturbations of large magnitude but small spatial extent in elliptic PDE problems and reaction-diffusion systems. This theory has been used to study mean first passage times in the context of molecular and cellular biology [30, 31, 32, 33, 34], pattern formation in biological development [35, 36], ecological dynamics [37], and diffusive signaling and chemoreception [38, 39, 40, 41]. Finally, since the domain geometries considered in this article are irregular (in the sense of containing curved boundaries), we numerically discretize our spatial differential operators on scattered node sets using Radial Basis Function (RBF)-based finite difference (RBF-FD) formulas. RBF-FD formulas have been used in a variety of applications and are widely known to provide both numerical stability and high orders of convergence [42, 43, 44, 45, 46, 47, 48, 49, 50, 51].

The rest of the paper is organized as follows. We first briefly summarize our main results in the subsection below. In section 2, we give a precise and dimensionless formulation of our stochastic PDE model. In section 3, we use asymptotic analysis to study the mean of this stochastic PDE. Section 4 applies our general asymptotic results to some examples. We describe our numerical methods and use them to verify our asymptotic results in section 5. We conclude by discussing the biological relevance of our model and describing future work.

1.1 Main results

Let $c(x,t)$ be the concentration of a neurotransmitter in the extracellular space $E\subset\mathbb{R}^{3}$ in some region of the brain. We assume that the neurotransmitter diffuses in the extracellular space,

[TABLE]

Furthermore, we assume that a set of neurons projects to this brain region, so that the domain $E$ contains $N\geq 1$ holes in its interior which represent the corresponding nerve varicosities, see Figure 1.

Each neuron fires action potentials and thus switches between a quiescent state and a firing state. We suppose that neuron firing is controlled by a continuous-time, reversible Markov jump process. We make no specific assumptions about firing correlations between neurons, and thus the neurons may fire synchronously, independently, or with some nontrivial correlations. When a neuron is not firing (quiescent), it absorbs neurotransmitter. When a neuron is firing, it releases neurotransmitter. Therefore, we impose the following stochastically switching boundary conditions at nerve varicosities,

[TABLE]

where $\partial_{\nu}$ is the outward unit normal derivative and $\partial V_{n}$ is the boundary of the $n$-th nerve varicosity. Here, ${\varphi}>0$ is some constant flux of neurotransmitter into the extracellular space when the neuron is firing. We assume that the neurotransmitter does not leave the region through the outer boundary, and so we impose a no flux boundary condition there,

[TABLE]

Summarizing, our model consists of the diffusion equation on a three-dimensional domain with interior holes that randomly switch between being sources and sinks (this model was first considered in [10] with spherical varicosities).

Let $l_{0}$ denote the characteristic length scale describing the size of each nerve varicosity, and let $l_{1}>0$ denote the characteristic distance between varicosities. Assuming $\varepsilon:=l_{0}/l_{1}\ll 1$, we derive the following asymptotic expansion for the large time expected neurotransmitter concentration,

[TABLE]

The $\mathcal{O}(\varepsilon)$ term is

[TABLE]

where $\pi_{1}^{(n)}\in(0,1)$ is the proportion of time that the $n$-th nerve varicosity is firing, $\pi_{0}^{(n)}=1-\pi_{1}^{(n)}$, and $C_{n}$ and $D_{n}$ are dimensionless constants determined by the shape of the $n$-th varicosity (see (27) and (36)). Notice that the $\mathcal{O}(\varepsilon)$ term is spatially constant and independent of the following factors: (i) the positions of the nerve varicosities, (ii) the number of nerve varicosities (provided the proportion of varicosities of a certain shape and firing probability is fixed), (iii) the size and geometry of the extracellular space, (iv) firing correlations between the neurons, and (v) the rate of switching between firing and quiescent states. To see point (ii), notice that $\theta^{(1)}$ can be written as

[TABLE]

where $\mu_{D}$ and $\mu_{C}$ are the averages,

[TABLE]

Hence, if the number $N$ of nerve varicosities changes, but the proportion of varicosities of a certain shape and firing probability is fixed, then $\mu_{D}$ and $\mu_{C}$ are unchanged, and thus $\theta^{(1)}$ is unchanged.

The $\mathcal{O}(\varepsilon^{2})$ term in (1) shows how these factors affect the average neurotransmitter concentration. This term is the sum of a constant $\theta^{(2)}$ and function $\Sigma(x)$. For simplicity, we omit the formulas for $\theta^{(2)}$ and $\Sigma(x)$ from this section (we give the formula in section 3 below, see (49)-(51)). However, we note that if all the varicosities have the same firing probability and shape ($\pi_{1}^{(n)}=\pi_{1}^{(m)}$, $C_{n}=C_{m}$, and $D_{n}=D_{m}$ for all $n,m\in\{1,\dots,N\})$, then $\Sigma(x)\equiv 0$, and thus this term is spatially constant.

2 Problem setup

We now describe a dimensionless version of the model in section 1.1 in further detail. Let $U\subset\mathbb{R}^{3}$ be open, connected, and bounded with a smooth boundary, $\partial U$. To represent the varicosities (or holes), let $\{V_{n}\}_{n=1}^{N}$ be $N$ open subsets of $U$, each with smooth boundary $\partial V_{n}$ and “radius” $\mathcal{O}(\varepsilon)$. That is, we assume that there exists $N$ points $\{x_{n}\}_{n=1}^{N}$ in $U$ such that $V_{n}\to x_{n}$ uniformly as $\varepsilon\to 0$. We assume that the holes are well-separated in the sense that

[TABLE]

and

[TABLE]

We define the extracellular space, $E$, as the domain minus the holes,

[TABLE]

To describe the state of each neuron (either firing or quiescent), let $\{J(t)\}_{t\geq 0}$ be an irreducible continuous-time Markov jump process on a finite state space $\mathcal{J}\subseteq\{0,1\}^{N}$, where $J_{n}(t)\in\{0,1\}$ denotes its $n$-th component for $n\in\{1,\dots,N\}$. That is, $J(t)$ is the $N$-dimensional vector,

[TABLE]

with $J_{n}(t)\in\{0,1\}$ for $n\in\{1,\dots,N\}$. We emphasize that the state space $\mathcal{J}$ depends on the firing correlations between neurons. For example, if the neurons fire synchronously, then $J$ consists of only two elements, namely the vector of length $N$ with all [math]’s and the vector of length $N$ with all $1$’s,

[TABLE]

On the other hand, if the neurons fire independently, then $\mathcal{J}=\{0,1\}^{N}$.

Let $Q\in\mathbb{R}^{|\mathcal{J}|\times|\mathcal{J}|}$ denote the infinitesimal generator matrix of $J(t)$. Letting $Q_{\mathbf{i},\mathbf{j}}\in\mathbb{R}$ denote the entry in the row corresponding to state $\mathbf{i}\in\mathcal{J}$ and the column corresponding to state $\mathbf{j}\in\mathcal{J}$, recall that $Q_{\mathbf{i},\mathbf{j}}\geq 0$ gives the rate that $J$ jumps from state $\mathbf{i}\in\mathcal{J}$ to state $\mathbf{j}\in\mathcal{J}$ for $\mathbf{i}\neq\mathbf{j}$, and the diagonal entries $Q_{\mathbf{i},\mathbf{i}}$ are chosen so that $Q$ has zero row sums [52].

Irreducibility ensures that $J(t)$ has a unique invariant distribution. That is, there exists a unique $\rho\in\mathbb{R}^{|\mathcal{J}|}$ satisfying $Q^{T}\rho=0$ and $\sum_{\mathbf{j}\in\mathcal{J}}\rho_{\mathbf{j}}=1$, where $Q^{T}$ denotes the transpose of $Q$. We further assume that $Q$ is diagonalizable with all real eigenvalues. We note that if $J(t)$ is reversible (satisfies detailed balance), then these assumptions on $Q$ are satisfied. To see this, note that reversibility means

[TABLE]

Hence, a straightforward calculation shows that $Q$ is self-adjoint with respect to the inner product

[TABLE]

The spectral theorem thus ensures that $Q$ is diagonalizable with all real eigenvalues.

We note that assuming $Q$ is diagonalizable with all real eigenvalues is slightly restrictive. For example, if $N=3$ neurons always fire in a prescribed order (the first neuron fires, then the second, then the third, then the first, etc.), then a quick calculation shows that $Q$ has complex eigenvalues. However, $Q$ is diagonalizable with all real eigenvalues in many cases. For example, if the neurons fire synchronously or independently, then $J(t)$ is reversible (and thus $Q$ is diagonalizable with all real eigenvalues by the argument above). In fact, if the $N$ neurons can be partitioned into subgroups so that each subgroup fires synchronously and different subgroups fire independently, then $J(t)$ is reversible.

Assume that $\{c(x,t)\}_{t\geq 0}$ is the $L^{2}(E)$-valued stochastic process that satisfies the diffusion equation in $E$ with a no flux condition at the outer boundary,

[TABLE]

and boundary conditions at each $V_{n}$ that switch according to the jump process $J_{n}(t)$,

[TABLE]

Since $c(x,t)$ is piecewise deterministic, it is straightforward to construct $c(x,t)$ by composing the deterministic solution operators corresponding to the $|\mathcal{J}|$-many sets of boundary conditions according to the path of $\{J(t)\}_{t\geq 0}$ (see [53] for a construction of a solution to a similar switching PDE).

This model is obtained from the model in section 1.1 by defining dimensionless time, length, and concentration variables, $\frac{D}{l_{1}^{2}}t$, $\frac{1}{l_{1}}x$, and $l_{1}^{3}c$, for some length scale, $l_{1}>0$. In these variables, the dimensionless inhomogeneous flux boundary condition (when a neuron is firing) is $l_{1}^{4}{\varphi}/D$, but we have taken this to be unity in (7) without loss of generality. In particular, if $c$ denotes the solution to (4)-(7), then the solution to the problem with $l_{1}^{4}{\varphi}/D$ on the righthand side of (7) is simply $(l_{1}^{4}{\varphi}/D)c$ by linearity. Hence, the dimensional result in (1) is obtained by solving (4)-(7), then multiplying by $l_{1}^{4}{\varphi}/D$, and then multiplying by $l_{1}^{-3}$ to obtain a dimensional concentration.

3 Asymptotic analysis of mean neurotransmitter

3.1 Decompose and diagonalize mean

To analyze the large time mean of $c(x,t)$, we decompose the mean based on the state of the jump process $J$. Specifically, for each element $\mathbf{j}=(j_{1},\dots,j_{N})\in\mathcal{J}\subseteq\{0,1\}^{|\mathcal{J}|}$ of the state space of $J$, define

[TABLE]

Here, $1_{J(t)=\mathbf{j}}$ denotes the indicator function,

[TABLE]

Putting the deterministic functions in (8) into a vector, $\mathbf{v}(x)=\{v_{\mathbf{j}}\}_{\mathbf{j}\in\mathcal{J}}\in\mathbb{R}^{|\mathcal{J}|}$, we have by Theorem 6 in [10] (see also [21]) that

[TABLE]

where $E:=U\backslash\overline{\cup_{n=1}^{N}V_{n}}$ is $U$ minus the holes. Since there is a no flux boundary condition at the outer boundary regardless of the state of $J(t)$, we have that

[TABLE]

One advantage of the decomposition in (8) is that $v_{\mathbf{j}}$ must satisfy the boundary condition at nerve varicosities corresponding to $J(t)=\mathbf{j}\in\mathcal{J}$,

[TABLE]

where $\rho=\{\rho_{\mathbf{j}}\}_{\mathbf{j}\in\mathcal{J}}\in\mathbb{R}^{|\mathcal{J}|}$ is the unique invariant probability measure of $J(t)$. The boundary conditions (10)-(11) follow from Theorem 6 in [10].

Writing (11) in matrix notation, we have that

[TABLE]

where $I\in\mathbb{R}^{|\mathcal{J}|\times|\mathcal{J}|}$ is the identity matrix and $B_{n}\in\{0,1\}^{|\mathcal{J}|\times|\mathcal{J}|}$ is the diagonal matrix with entries,

[TABLE]

That is, $B_{n}$ is the $|\mathcal{J}|\times|\mathcal{J}|$ diagonal matrix whose diagonal entry corresponding to state $\mathbf{j}\in\mathcal{J}$ is 1 (respectively 0) if the $n$-th neuron is firing (respectively quiescent) when $J(t)=\mathbf{j}$. We emphasize that (12) is merely a compact way of writing the homogeneous Dirichlet conditions and inhomogeneous Neumann conditions that are given in (11). Hence, while (12) may look like a Robin condition, it is in fact a vector of homogeneous Dirichlet conditions and inhomogeneous Neumann conditions due to the structure of $B_{n}$ given in (13).

Notice that $\mathbf{v}$ satisfies the coupled PDE (9) with decoupled boundary conditions (10) and (12). We now decouple the PDE at the cost of coupling the boundary conditions in (12). By assumption, we can diagonalize $Q^{T}$,

[TABLE]

where

[TABLE]

and the first column of $P\in\mathbb{R}^{|\mathcal{J}|\times|\mathcal{J}|}$ is the invariant measure $\rho\in\mathbb{R}^{|\mathcal{J}|}$, and the first row of $P^{-1}\in\mathbb{R}^{|\mathcal{J}|\times|\mathcal{J}|}$ is all 1’s. Hence, setting

[TABLE]

we have that $\mathbf{u}$ satisfies

[TABLE]

with no flux boundary conditions at the outer boundary

[TABLE]

and boundary conditions at the interior holes given by

[TABLE]

3.2 Inner and outer expansions

Using (14), observe that the first component of $\mathbf{u}$ (denoted by $(\mathbf{u})_{1}$) is the large time mean

[TABLE]

Therefore, it remains to solve (15) subject to (16) and (17) in order to find the large time mean. We approximate the solution to this ${|\mathcal{J}|}$-dimensional boundary value problem in the case that the nerve varicosities are small, $\varepsilon\ll 1$. Specifically, we construct inner (or local) solutions that are valid in an $\mathcal{O}(\varepsilon)$ neighborhood of each varicosity and then match these to an outer (or global) solution that is valid away from these neighborhoods, see Figure 2. The asymptotic analysis that follows is formal and is related to the strong localized perturbation analysis pioneered in [27].

The outer problem is

[TABLE]

Notice that from the perspective of the outer solution, the $N$ holes (varicosities), $\{V_{n}\}_{n=1}^{N}$, have shrunk to the points, $\{x_{n}\}_{n=1}^{N}$. To construct an inner solution that is valid near the $n$-th nerve varicosity, we introduce the stretched coordinate

[TABLE]

and define

[TABLE]

This inner solution satisfies

[TABLE]

where $U_{n}:=\varepsilon^{-1}V_{n}$ is the magnified hole (note that $y\in\partial U_{n}$ if and only if $x\in\partial V_{n})$. Due to the structure of $B_{n}$ in (13), note that the factor of $\varepsilon$ in the first term in the boundary condition in (19) does not affect the solution. Hence, in the rest of the calculation we replace (19) by

[TABLE]

since this does not affect the solution $\mathbf{w}_{n}$.

Next, we expand these inner and outer solutions in powers of $\varepsilon$,

[TABLE]

It follows that the $k$-th term in the outer solution expansion satisfies

[TABLE]

Similarly, the $k$-th term in the inner solution expansion satisfies

[TABLE]

where $\mathbf{w}_{n}^{(-k)}\equiv 0$ if $k>0$ and $1_{k=1}$ is 1 if $k=1$ and [math] otherwise. The matching condition between the outer and inner solutions is that the near-field behavior of the outer solution as $x\to x_{n}$ must agree with the far-field behavior of the inner solution as $|y|\to\infty$. We write this formally as

[TABLE]

3.3 $\mathcal{O}(\varepsilon^{0})$ term

Consider the problem with no holes in the domain, $\varepsilon=0$. Since $\lambda_{1}=0$ (the first diagonal entry of $\Lambda$), the unperturbed outer solution, $\mathbf{u}^{(0)}(x)$, may be nontrivial in the first component, and we have the form

[TABLE]

for some constant $\theta^{(0)}\in\mathbb{R}$ to be determined.

With this information about $\mathbf{u}^{(0)}$, we can solve the inner problems for $\mathbf{w}_{n}^{(0)}$, $n\in\{1,\dots,N\}$. Let $f(y)$ be the real-valued solution to

[TABLE]

Then, the solution to the inner problem near the $n$-th hole is

[TABLE]

To see this, observe that (26) satisfies (22) since $f$ is harmonic. Further, (26) satisfies (23) since by the structure of $B_{n}$ in (13),

[TABLE]

and $f=1$ on $\partial U_{n}$. Finally, (26) satisfies (24) since $f\to 0$ as $|y|\to\infty$.

We can now use (26) to determine the leading order behavior of $\mathbf{u}^{(1)}$ as $x\to x_{n}$. The function $f(y)$ can be found in closed form only for a few simples shapes, but for a general shape it has the far-field asymptotic behavior

[TABLE]

where $C_{n}>0$ is called the electrostatic capacitance of ${U}_{n}$, and $\Phi_{n}\in\mathbb{R}^{3}$ is its dipole vector, and both are determined by $U_{n}$ [54]. In general, $C_{n}$ must be calculated numerically, but it is known analytically for certain shapes, see Table 1 in [30]. In particular, if $U_{n}$ is a sphere, then $C_{n}$ is its radius. We will not need the value of $\Phi_{n}$.

Hence, (26) and (27) imply that $\mathbf{w}_{n}^{(0)}$ has the far-field behavior

[TABLE]

where

[TABLE]

Plugging (28) into the matching condition (24) and using $y=\varepsilon^{-1}|x-x_{n}|$, we have that

[TABLE]

Hence, the outer solution $\mathbf{u}^{(1)}$ has the singular behavior as $x\to x_{n}$,

[TABLE]

We emphasize that $\mathbf{u}^{(1)}(x)$ is part of the expansion of the outer solution and is thus valid for $x$ outside an $\mathcal{O}(\varepsilon)$ neighborhood of each of the holes. This singular behavior can be incorporated into the problem (20) for $\mathbf{u}^{(1)}$ in terms of the Dirac delta function as

[TABLE]

Integrating the first component of the $|\mathcal{J}|$-dimensional vector PDE in (30) and using the divergence theorem and the fact that the first row of $\Lambda$ is all zeros, we obtain that the sum over $n$ of the first components of $\chi_{n}^{(1)}$ must be zero,

[TABLE]

Now, using (13), (14), (25), and (29), we obtain

[TABLE]

Therefore,

[TABLE]

and thus $\mathbf{w}_{n}^{(0)}=0$ by (26).

3.4 $\mathcal{O}(\varepsilon^{1})$ term

Since $\theta^{(0)}=0$, we have that $\chi_{n}^{(1)}=0$ for all $n\in\{1,\dots,N\}$ from (29). Thus, the solution to (30) is the vector

[TABLE]

for some constant $\theta^{(1)}\in\mathbb{R}$ to be determined.

With this information about $\mathbf{u}^{(1)}$, we can solve the inner problem for $\mathbf{w}_{n}^{(1)}$, $n\in\{1,\dots,N\}$. Let $g(y)$ be the real-valued solution to

[TABLE]

The solution to the inner problem near the $n$-th hole is

[TABLE]

To see this, observe that (34) satisfies (22) since $f$ and $g$ are harmonic. Further, (34) satisfies (23) since

[TABLE]

and $f=\partial_{\nu}g=1$ on $\partial U_{n}$. Finally, (34) satisfies (24) since $f\to 0$ and $g\to 0$ as $|y|\to\infty$.

We can now use (34) to determine the leading order behavior of $\mathbf{u}^{(2)}$ as $x\to x_{n}$. Using the known far-field behavior in (27) for solution to the analogous Dirichlet problem, we assume that the function $g(y)$ has the far-field asymptotic behavior

[TABLE]

where $D_{n}>0$ and $\Psi_{n}\in\mathbb{R}^{3}$ are determined by the shape of $U_{n}$. In fact, integrating the PDE in (33) over a large sphere containing $U_{n}$ and using that $\partial_{\nu}g=1$ on $\partial U_{n}$ and (35), we obtain

[TABLE]

where $|\partial U_{n}|$ is the surface area of the magnified hole, $U_{n}:=\varepsilon^{-1}V_{n}$. We note that if $U_{n}$ is a sphere, then it thus follows that $D_{n}=C_{n}^{2}$, where $C_{n}$ is the radius of $U_{n}$. We will not need the value of $\Psi_{n}$.

Hence, $\mathbf{w}_{n}^{(1)}$ has the far-field behavior

[TABLE]

where

[TABLE]

Plugging (37) into the matching condition (24) and using $y=\varepsilon^{-1}|x-x_{n}|$, we have that

[TABLE]

Hence, $\mathbf{u}^{(2)}$ has the singular behavior as $x\to x_{n}$,

[TABLE]

This singular behavior can be incorporated into the problem (20) for $\mathbf{u}^{(2)}$ in terms of the Dirac delta function as

[TABLE]

Integrating the first component of the $|\mathcal{J}|$-dimensional vector PDE in (40) and using the divergence theorem and the fact that the first row of $\Lambda$ is all zeros, we obtain

[TABLE]

Now, using (13), (14), (32), and (38), we have that

[TABLE]

where we have defined the marginal probability that $J_{n}$ is 1 or 0,

[TABLE]

Therefore,

[TABLE]

Equation (42) generalizes a result proven in [10] using probabilistic methods for the special case that each hole is sphere.

3.5 $\mathcal{O}(\varepsilon^{2})$ term

The solution to (40) can be written in terms of Green’s functions. Let $G(x,\xi)$ denote the Neumann Green’s function, which satisfies

[TABLE]

where $|U|$ denotes the volume of the domain $U$. For $\lambda>0$, let $G_{\textup{H}}(x,\xi;\lambda)$ denote the modified Helmholtz Green’s function, which satisfies

[TABLE]

The functions $R$ and $R_{\textup{H}}$ are called the regular parts of the Green’s functions. In the special case that the domain $U$ is a sphere, $G$ and $G_{\textup{H}}$ are known explicitly [55, 30]. We take advantage of these explicit formulas in see section 4 below.

Putting these Green’s functions into a matrix, let $\mathbf{G}(x,\xi)$ denote the $|\mathcal{J}|\times|\mathcal{J}|$ diagonal matrix,

[TABLE]

where $(0,\lambda_{2},\dots,\lambda_{|\mathcal{J}|})$ are the diagonal entries of $\Lambda$. Hence, the solution to (40) is

[TABLE]

where $\Theta^{(2)}\in\mathbb{R}^{|\mathcal{J}|}$ is a vector of all zeros, except for the first component which is some constant $\theta^{(2)}\in\mathbb{R}$ to be determined.

Next, we use (45) to expand $\mathbf{u}^{(2)}$ as $x\to x_{n}$ and obtain

[TABLE]

where

[TABLE]

where $\Lambda$ is the diagonal matrix in (14) and where we have defined the $|\mathcal{J}|\times|\mathcal{J}|$ diagonal matrix,

[TABLE]

containing the regular parts of the Green’s functions.

Having obtained this asymptotic behavior of $\mathbf{u}^{(2)}$, we now return to the matching condition (24) to determine the far-field behavior of the inner solution, $\mathbf{w}_{n}^{(2)}$, $n\in\{1,\dots,N\}$. In particular, we find that

[TABLE]

Hence, since $\mathbf{w}_{n}^{(2)}$ must also satisfy (22)-(23), we have that

[TABLE]

Therefore, we obtain the next term in the far-field behavior of $\mathbf{w}_{n}^{(2)}$,

[TABLE]

where

[TABLE]

Plugging this into the matching condition (24) (similar to the procedure that led to (39)), we find that $\mathbf{u}^{(3)}$ has the leading-order singular behavior as $x\to x_{n}$,

[TABLE]

This singular behavior can be incorporated into the problem (20) for $\mathbf{u}^{(3)}$ in terms of the Dirac delta function as

[TABLE]

with no flux boundary conditions, $\partial_{\nu}\mathbf{u}^{(3)}=0,\;x\in\partial U$. Integrating the first component of the $|\mathcal{J}|$-dimensional vector PDE in (48) and using the divergence theorem and the fact that the first row of $\Lambda$ is all zeros, we obtain

[TABLE]

Using (13), (14) and (47), we solve this equation for $\theta^{(2)}$ and find

[TABLE]

3.6 Summary

We now put these calculations together and make some observations. In the small hole limit, $\varepsilon\to 0$, the large time expected neurotransmitter concentration is given asymptotically in the outer region $|x-x_{n}|\gg\mathcal{O}(\varepsilon)$ for $n\in\{1,\dots,N\}$ by

[TABLE]

where $\theta^{(1)}$ and $\theta^{(2)}$ are given by (42) and (49),

[TABLE]

\big{(}\chi_{n}^{(2)}\big{)}_{1} is (41), and $G(x,x_{n})$ is the Neumann Green’s function (43).

From (50), first observe that the large time expected neurotransmitter concentration is constant in space at order $\varepsilon$. Next, observe from (42) that this constant is independent of (i) the positions of the nerve varicosities, $\{x_{n}\}_{n=1}^{N}$, (ii) number $N$ of nerve varicosities (provided the proportion of varicosities of a certain shape and firing probability is fixed, see (2)), (iii) the size and geometry of the extracellular space, $E$, and (iv) firing correlations between the neurons. Elaborating on (iv), this means that the constant is unchanged if the neurons fire synchronously (perfectly correlated), independently, or with other correlations. Finally, (42) shows that this constant depends on the proportion of time spent firing but not on the overall rate of switching between firing and quiescent states.

From (50) and (41), further observe that if all the holes have a common shape (up to rotations) and firing probability, then \big{(}\chi_{n}^{(2)}\big{)}_{1}=0. Therefore in this case, observe from (50) that the large time expected neurotransmitter concentration is in fact constant in space at order $\varepsilon^{2}$. Furthermore, from (50), (49), and (46), observe that the order $\varepsilon^{2}$ term describes how the number, positions, firing rates, and correlations of the neurons and the geometry of the extracellular space affects the mean neurotransmitter concentration.

4 Examples

In this section, we consider some specific examples of firing correlations and geometries in order to illustrate and verify the general results of section 3.

4.1 Synchronous firing

First suppose that all $N$ neurons fire synchronously. That is, suppose $J(t)\in\mathcal{J}$ where the state space $\mathcal{J}$ consists of the two elements (i.e. $|\mathcal{J}|=2$) given in (3). Namely, $\mathcal{J}$ consists of the $N$-dimensional vector of all 0’s corresponding to every neuron being queiscent and the $N$-dimensional vector of all 1’s corresponding to every neuron firing. If we label these two states by [math] and $1$, then $J(t)$ has transition rates

[TABLE]

In this case, (9)-(11) simplify to the pair of PDEs,

[TABLE]

with boundary conditions

[TABLE]

It follows from (52) that

[TABLE]

If we further suppose that each $V_{n}$ is a sphere of radius $\varepsilon>0$ so that $C_{n}=D_{n}=1$ for $n\in\{1,\dots,N\}$, then

[TABLE]

and thus

[TABLE]

It then follows from (49) that

[TABLE]

Notice that our expression for $\theta^{(2)}$ (and thus our expression for the large time mean neurotransmitter in (50)) involves the Neumann Green’s function of the modified Helmholtz equation (44). Since this Green’s function cannot be computed analytically for a general domain $U\subset\mathbb{R}^{3}$, we now suppose that the domain $U$ is a sphere of radius $R_{0}>0$. In this case, the modified Helmholtz Green’s function satisfying (44) is given explicitly by [55]

[TABLE]

where

[TABLE]

Here, $I_{\nu}(x)$ (respectively $K_{\nu}(x)$) is the modified Bessel function of the first (respectively second) kind of order $\nu$, $p_{n}(\cos\gamma)$ is the Legendre polynomial of degree $n$ with argument $\cos\gamma=\cos\theta\cos\theta_{0}+\sin\theta\sin\theta_{0}\cos(\varphi-\varphi_{0})$, and $(r,\theta,\varphi)$ (respectively $(r_{0},\theta_{0},\varphi_{0})$) denotes the spherical coordinates of $x$ (respectively $x_{0}$).

Summarizing, we have obtained the formula

[TABLE]

where $\theta^{(1)}$ and $\theta^{(2)}$ have the analytical formulas (56) and (57), where the Green’s functions are given in (58)-(59). In section 5 below, we show that (60) agrees with an independent numerical approximation to (53)-(54).

4.2 Spherically symmetric

We now consider a simple example, which is the only example in which we can find an analytical formula for the mean neurotransmitter. We then use this example to further check our asymptotic result in (60).

Suppose there is a single spherical hole of radius $\varepsilon>0$ located in the center of a sphere of radius $R_{0}>\varepsilon$. That is, $N=1$, with

[TABLE]

Let $\alpha$ and $\beta$ be as in the previous example. Finding the mean neurotransmitter reduces to solving (53) in the annulus $|x|\in(\varepsilon,R_{0})$ with boundary conditions

[TABLE]

It is straightforward to use radial symmetry to find the constant value of the large-time mean neurotransmitter. In particular, letting $r=|x|$, diagonalizing the problem as in (15)-(17), using the matrices in (55), and using radial symmetry gives that

[TABLE]

satisfies

[TABLE]

with boundary conditions

[TABLE]

It follows directly from (62) and the first boundary condition in (64) that $u_{0}$ is spatially constant. We can then solve (63), the second boundary condition in (64), and the second boundary condition in (65) to obtain $u_{1}$ (since $\tfrac{d}{dr}u_{0}=0$). From this value of the function $u_{1}$, we then obtain the constant value of $u_{0}$ from the first boundary condition in (65), which yields

[TABLE]

where

[TABLE]

With this explicit expression, one can obtain the Taylor series expansion for $\varepsilon\ll 1$,

[TABLE]

We now verify that our more general asymptotic result in (60) reproduces (67) in this special case. It is immediate that $\theta^{(1)}$ in (56) agrees with (67). Furthermore, in this special case we have that $\theta^{(2)}$ in (57) reduces to

[TABLE]

Now, it follows from (59) that

[TABLE]

To derive (69), let $r=|x_{1}|$ and note that (59) yields

[TABLE]

Now, it is a property of the modified Bessel function that

[TABLE]

Hence, taking the limit $|x_{1}|=r\to 0+$ in (70) yields only the $n=0$ term,

[TABLE]

Plugging the value $a_{0}=K_{3/2}(\sqrt{\alpha+\beta}R_{0})/I_{3/2}(\sqrt{\alpha+\beta}R_{0})$ from (59) into (71) and simplifying yields (69).

Plugging (69) into (68) yields that

[TABLE]

which agrees with (67).

4.3 Independent firing

Suppose that each of the $N\geq 1$ nerve varicosities, ${V}_{n}$, fires according to an independent Markov jump process $J_{n}(t)$ with transition rates $\alpha_{n},\beta_{n}$,

[TABLE]

and $J(t)=(J_{1}(t),\dots,J_{N}(t))\in\mathcal{J}=\{0,1\}^{N}$. In this case,

[TABLE]

where $\delta_{i,j}$ denotes the Kronecker delta and $q_{i,j}^{(n)}$ are the entries in the generator $Q^{(n)}$ of the Markov process $J_{n}(t)$ controlling the firing of the $n$-th nerve varicosity,

[TABLE]

Observe that $Q_{\mathbf{i}\mathbf{j}}=0$ if $\mathbf{i}$ and $\mathbf{j}$ differ in more than one component, and either $Q_{\mathbf{i}\mathbf{j}}=\alpha_{n}$ or $Q_{\mathbf{i}\mathbf{j}}=\beta_{n}$ if $\mathbf{i}$ and $\mathbf{j}$ differ in exactly the $n$-th position. Further, $Q_{\mathbf{i}\mathbf{i}}=\sum_{n=1}^{N}q_{i_{n},i_{n}}<0$ with $q_{i_{n},i_{n}}=-(1-i_{n})\beta_{n}-i_{n}\alpha_{n}$. Also,

[TABLE]

Since the components of $J(t)$ are independent, the probability that the $n$-th component is 1 is $\frac{\alpha_{n}}{\alpha_{n}+\beta_{n}}$ and thus

[TABLE]

If $\alpha_{m}=\alpha_{n}$ and $\beta_{m}=\beta_{n}$ for all $m,n\in\{1,\dots,N\}$ then $\theta^{(1)}$ is the same as if the neurons fired synchronously (as noted above in section 3.6).

To simplify the calculation of $\theta^{(2)}$, we assume that there are only two varicosities, $N=2$, and that each one is a sphere of radius $\varepsilon>0$ so that $C_{n}=D_{n}=1$ for $n=1,2$. Further suppose that each varicosity has the same firing rate so that $\alpha_{n}=\alpha$ and $\beta_{n}=\beta$ for $n=1,2$. In this case, $\theta^{(1)}=\beta/\alpha$. Furthermore, it follows that

[TABLE]

and

[TABLE]

Using (49) and (72), a tedious but straightforward calculation yields

[TABLE]

where $x_{1}$ and $x_{2}$ are the centers of the two varicosities.

We now compare the case of independent firing considered here to the case of synchronous firing in section 4.1 above. Let $\theta_{\text{indep}}^{(2)}$ denote the value of $\theta^{(2)}$ in (73) for independent firing and $\theta_{\text{sync}}^{(2)}$ denote the value of $\theta^{(2)}$ in (57) for synchronous firing when $N=2$. Comparing (57) to (73), we have that

[TABLE]

To understand (74), observe that the Green’s function $G_{\textup{H}}(x,\xi;\lambda)$ for $\lambda>0$ is the steady-state concentration at position $x\in U$ of a diffusing chemical that degrades at rate $\lambda>0$, has a source at position $\xi\in U$, and reflects from the boundary $\partial U$. It follows that (i) $G_{\textup{H}}(x,\xi;\lambda)>0$, (ii) $G_{\textup{H}}(x,\xi;\lambda)$ grows as $x$ approaches $\xi$, and (iii) $G_{\textup{H}}(x,\xi;\lambda)$ grows as $\lambda$ decreases (these three intuitive properties can be proved by analyzing (44)).

Using these three properties of $G$, equation (74) implies the following three conclusions. First, the neurotransmitter concentration for synchronous firing is strictly greater than the neurotransmitter concentration for independent firing. Second, this difference between synchronous and independent firing is more pronounced when the varicosities are close to each other. Third, this difference is also more pronounced when the overall rate of firing, $\alpha+\beta$, is slow. Moreover, while (74) concerns perfectly synchronous versus perfectly independent firing for $N=2$ varicosities, we conjecture that these three conclusions hold when comparing more correlated versus less correlated firing of $N\geq 2$ varicosities.

5 Numerical analysis of mean neurotransmitter

In this section, we verify our asymptotic results through numerical simulations. For these simulations, we consider the case of synchronous firing described in sections 4.1-4.2 above for three different domain geometries. First, we consider the spherically symmetric case of section 4.2. Since the large-time mean neurotransmitter can be found explicitly in this case (see (66)), this case serves to validate our numerical methods. We then move on to more complicated geometries that cannot be solved analytically to verify our asymptotic result in (60). Before detailing these numerical experiments and their results, we first describe the discretization of the PDE (53) and boundary conditions (54) on these domains.

5.1 Spatial discretization

As mentioned in section 1, we use the RBF-FD method for spatial discretization. More specifically, we use a specific variant of the RBF-FD method (developed by the second author) called the overlapped RBF-FD method [48, 49]. This variant of the RBF-FD method uses a judicious mix of centered and one-sided differences on scattered nodes, leading to fewer finite difference stencils than collocation points. This allows for improved computational efficiency over standard RBF-FD without sacrificing high-order accuracy or numerical stability. We now briefly describe the overlapped RBF-FD method for approximating the Laplacian; the same discussion applies to any linear operator (including the Neumann boundary condition operator).

Let $X=\{{{x}}_{k}\}_{k=1}^{N}$ be a set of collocation nodes on some domain $\Omega\subset\mathbb{R}^{3}$ with boundary $\partial\Omega$. Define the stencil $P_{k}$ to be the set of nodes containing the nodes ${{x}}_{{\mathcal{I}}^{k}_{1}}$, $\{{{x}}_{{\mathcal{I}}^{k}_{2}},\ldots,{{x}}_{{\mathcal{I}}^{k}_{n}}\}$, i.e., ${{x}}_{{\mathcal{I}}^{k}_{1}}$ and its $n-1$ nearest neighbors. Here, ${\mathcal{I}}^{k}_{j}$, $j=1,\ldots n$ are indices that map stencil nodes into the set $X$. Without loss of generality, we confine our discussion to the stencil $P_{1}$ and the Laplacian $\Delta$. First, given a parameter $\delta\in(0,1]$, we can define a ball associated with stencil $P_{1}$ with center ${{x}}_{{\mathcal{I}}^{1}_{1}}$ and radius $r_{1}$ given by:

[TABLE]

This allows us to define the subset $\mathbb{B}_{1}\subset P_{1}$ as:

[TABLE]

i.e., $\mathbb{B}_{1}$ contains all the nodes in $P_{1}$ that lie within the ball of radius $r_{1}$ centered at ${{x}}_{{\mathcal{I}}^{1}_{1}}$. Finally, let the global indices of the $p_{1}$ nodes in $\mathbb{B}_{1}$ be given by:

[TABLE]

The overlapped RBF-FD method now involves using all the nodes in $P_{1}$ (${{x}}\in P_{1}$) to compute RBF-FD weights for the nodes in $\mathbb{B}_{1}$ (${{y}}\in\mathbb{B}_{1}$). To do so, we first define the augmented local RBF interpolant on $P_{1}$:

[TABLE]

where $\|{{x}}-{{x}}_{{\mathcal{I}}^{1}_{j}}\|^{m}$ is the polyharmonic spline (PHS) RBF of degree $m$ ($m$ is odd), and $\psi^{1}_{i}({{x}})$ are the $M$ monomials corresponding to the total degree polynomial of degree $\ell$ in $d$ dimensions. In order to compute the RBF-FD weights $w^{1}_{j}({{y}})$ for the Laplacian, we impose the following two sets of conditions:

[TABLE]

The interpolant (78) and the two conditions (79)–(80) can be collected into the following block linear system:

[TABLE]

where

[TABLE]

$W_{1}$ is the matrix of overlapped RBF-FD weights, with each column containing the RBF-FD weights for a point ${{y}}\in\mathbb{B}_{1}$. The linear system (81) has a unique solution if the nodes in $P_{1}$ are distinct [56, 57]. The matrix of polynomial coefficients $\Lambda_{1}$ is a set of Lagrange multipliers that enforces the polynomial reproduction constraint (80), and can be discarded. This procedure is repeated for all points ${{x}}\in X$ that were not in the sets $\mathbb{B}$. The weights are then collected into the rows of a sparse matrix $L$, with each stencil populating multiple rows of $L$. For more details, see [48, 49]. In this setting, for a given function $f$, the error in approximating $\Delta f$ at a point ${{y}}\in\mathbb{B}$ is given by [58]:

[TABLE]

where $h({{y}})$ is the largest distance between the points ${{y}}$ and all points in a given stencil; clearly, the error depends on the polynomial degree $\ell$. Thus, if we wished for a convergence rate of order $\xi$, we would need to set the polynomial degree as $\ell=\xi+1$. Following [49], we set the stencil size $n$ to be $n=2M+\left\lfloor\ln(2M)\right\rfloor$, where $M={\ell+d\choose d}$. Finally, again from [49], we select $\delta$ to be:

[TABLE]

5.2 Numerical solution of the PDE

Having discussed the discretization of differential operators using overlapped RBF-FD, we now turn to the discretization of the PDE (53) and boundary conditions (54). In terms of the solution $v_{0},v_{1}$ to the PDE (53), recall that (18) implies that the large time expected neurotransmitter is

[TABLE]

We denote the numerical approximations of $v_{0}$, $v_{1}$, and $v$ by ${{\mathbb{V}_{0}}}$, ${{\mathbb{V}_{1}}}$, and ${{\mathbb{V}}}:={{\mathbb{V}_{0}}}+{{\mathbb{V}_{1}}}$.

Now, the problem (53)-(54) is a pair of coupled Helmholtz-type PDEs, with the variables sharing the same Neumann boundary condition on the outer domain boundary, and having different boundary conditions on any inner domain boundaries. Let $X=\{{{x}}_{k}\}_{k=1}^{N}$ be some set of nodes in the domain of interest; in our case, we used the node generator from [59] to generate these node sets. For convenience, we partition $X$ into two sets:

[TABLE]

where $X_{i}$ are the $N_{i}$ nodes in the interior of the domain, $X_{b}=[X_{b}^{o}\ X_{b}^{i}]$ are the $N_{b}$ boundary nodes; the superscripts $o$ and $i$ indicate outer and inner boundaries, respectively. To solve the PDE, we enforce that the PDEs hold up to and including the domain boundaries, and that the boundary conditions hold true at the boundaries. This allows us to use a set of ghost nodes outside the domain to enforce boundary conditions, giving us additional numerical stability [49]. Let $X_{g}=[X_{g}^{o}\ X_{g}^{i}]$ be the set of ghost nodes associated with the boundary nodes. Thus, in practice, we are operating on an extended node set $\tilde{X}$ given by

[TABLE]

Given this partitioning of the node sets, There are many reasonable possible orderings of the unknowns. For simplicity, we choose to order the unknowns as

[TABLE]

Note that

[TABLE]

that is, we order ${{\mathbb{V}_{0}}}$ at the boundary nodes as ${{\mathbb{V}_{0}}}$ at the outer boundary nodes followed by ${{\mathbb{V}_{0}}}$ inner boundary nodes; the same ordering is used for ${{\mathbb{V}_{1}}}$. Of course, for consistency, the same ordering is followed for ghost nodes as well. Now, given an arbitrary function $f$, we know that $L_{\tilde{X}}f_{\tilde{X}}\approx\left.\Delta f\right|_{\tilde{X}}$. Consistent with the above ordering, we partition the sparse matrix $L_{\tilde{X}}$ as follows:

[TABLE]

where the subscript $ii$ indicates mapping interior nodes to interior nodes, and so forth. We denote the discretized Neumann operator corresponding to the outer boundary by $D^{o}$ (a sparse matrix). We also require a matrix for Neumann boundary conditions on the inner boundaries; denote this by $D^{i}$. The Dirichlet boundary conditions on the interior do not require a special boundary operator, and can be enforced using identity matrices and zero matrices. We partition these matrices similar to $L$, with the exception that there are no rows mapping any nodes to interior points. To aid in the discretization of the PDE, define the boundary operators for ${{\mathbb{V}_{0}}}$ to be

[TABLE]

which corresponds to enforcing a Neumann boundary condition on the outer boundary, and a Dirichlet condition on the inner boundary. Similarly, the boundary operators for ${{\mathbb{V}_{1}}}$ are:

[TABLE]

which corresponds to Neumann boundary conditions on both outer and inner boundaries. With this notation in place, the fully discretized (block) system of PDEs and boundary conditions is given by:

[TABLE]

where $\tilde{\Theta}=[0\ \Theta]^{T}$, and $\Theta$ is a vector with length equal to the number of inner boundary points with each entry as the constant $\frac{\beta}{\alpha+\beta}$. In the above system, $I$ refers to the identity matrix, with the subscripts indicating the size. [math] in each case refers to a matrix of zeros; since the size of these zero blocks can be inferred from the rows in which they sit, we leave out subscripts for brevity. The block matrix $\tilde{L}$ is very sparse, since its individual blocks are themselves sparse matrices. Note that this matrix has dimensions $2(N+N_{b})$. Since the continuous problem is well-posed, a good discretization of the continuous problem should result in a matrix $\tilde{L}$ that does not have a nullspace. In practice, this translates into a requirement on both the node sets used to discretize the domain, and on the RBF-FD discretizations of the Laplacian and Neumann operators. We are able to invert the matrix $\tilde{L}$ using Matlab’s built-in sparse direct solver without any difficulties. We choose the sparse direct solver for simplicity, rather than using a preconditioned iterative solve. We leave such extensions for future work.

5.3 Numerical Results

In this section, we present two types of numerical results. First, we show results that validate our numerical method on the spherically symmetric problem of section 4.2 for which we have the analytical solution (66). Then, using this validation test to guide our choices, we present a numerical investigation of our asymptotic results in two different domain geometries.

5.3.1 Validation of numerical methods

We first consider the spherically symmetric problem of section 4.2 in which the domain is the unit sphere with a spherical hole of radius $\varepsilon>0$ centered at the origin. In this scenario, the large-time expected neurotransmitter $v=v_{0}+v_{1}$ is exactly constant in space and this constant is given explicitly in (66). To validate our numerical methods, we thus want to show that (i) the numerical solution is approximately constant in space and (ii) this approximately constant value of the numerical solution converges to the value of (66).

Before checking these two points, we note that since this is a problem in three space dimensions and spatial refinement is relatively expensive, we instead measure convergence by order refinement. In the overlapped RBF-FD method, this is accomplished by increasing $\xi$; we use $\xi\in[2,5]$ for this study. In this case, we expect the error to decrease as $\xi$ increases, regardless of the hole radius $\varepsilon$. However, when we first attempted this study for holes of decreasing radius, we observed that the errors were higher for smaller holes, primarily due to our use of a fixed $h$ in the spatial discretization. After some experimentation, we found that maintaining a fixed ratio of $h/\varepsilon$ was critical to ensuring that the spatial errors were independent of hole radius (we take $h/\varepsilon\approx 1/7$). We continue to use this refinement strategy in section 5.3.2.

To show that the numerical solution is approximately constant in space, the left panel of Figure 3 plots the following measure of spatial variation,

[TABLE]

where $\mathbb{V}=\mathbb{V}_{0}+\mathbb{V}_{1}$ is the numerical approximation of $v=v_{0}+v_{1}$ and the maximums and minimums are taken over all interior nodes. This plot confirms that the deviation from a spatially constant solution decays under order refinement. Indeed, this plot shows that the maximum and minimum of the solution are well within $1\%$ for the lowest order method and well within $0.1\%$ for the highest order method, and this holds regardless of hole radius $\varepsilon$.

To show that the (approximately constant) value of the numerical solution approaches the analytical value $v=v_{0}+v_{1}$ under order refinement, the right panel of Figure 3 plots the normalized error, $|v-\overline{\mathbb{V}}|/v$, where $\overline{\mathbb{V}}$ denotes the value of the numerical solution averaged over all interior grid points. This plot shows that the error drops off under order refinement, reaching much less than $0.1\%$ for the highest order method.

5.3.2 Investigation of asymptotics

Having validated our numerical methods in section 5.3.1, we employ our numerical discretization to verify some of our asymptotic results. Specifically, we verify the asymptotic result in (60) in section 4.1 in the following two geometries. The first is for a spherical domain of radius $R_{0}=1$ with two holes of radius $\varepsilon$ centered at the Cartesian coordinates $(-0.3,-0.3,-0.3)\in\mathbb{R}^{3}$ and $(0.3,0.3,0.3)\in\mathbb{R}^{3}$. The second is a spherical domain of radius $R_{0}=1.2$ with the same holes. To ensure that the discretization error is small, we use $\xi=4$ (though similar results were obtained with $\xi=3$ and $\xi=5$).

Now, the asymptotic result in (60) predicts that the solution is approximately constant outside an $\mathcal{O}(\varepsilon)$ neighborhood of each of the holes. Letting ${{\mathbb{V}}}={{\mathbb{V}_{0}}}+{{\mathbb{V}_{1}}}$ denote the numerical estimate to $v=v_{0}+v_{1}$ where $v_{0},v_{1}$ are the (unknown) exact solutions to (53)-(54), the left panel of Figure 4 confirms that the numerical solution is approximately spatially constant and that the deviation from constant decays as $\varepsilon$ decays. Indeed, the left panel plots the difference, $\max\mathbb{V}-\min\mathbb{V}$, where the maximum and minimum are taken over all interior grid points, and this difference is approximately $\mathcal{O}(\varepsilon^{3})$ or $\mathcal{O}(\varepsilon^{4})$, which agrees with the asymptotic result in (60).

The right panel of Figure 4 plots the absolute value of the difference between the asymptotic estimate, $\varepsilon\theta^{(1)}+\varepsilon^{2}\theta^{(2)}$, in (60) and the numerical solution averaged over all interior grid points, $\overline{\mathbb{V}}$. In agreement with (60), this difference is $\mathcal{O}(\varepsilon^{3})$. Indeed, the lines of best fit to the $R_{0}=1$ and $R_{0}=1.2$ curves in the semi-log plot of the right panel yields respective slopes of 3.19 and 2.96.

6 Discussion

We have used asymptotic analysis to study the large-time mean of a stochastic PDE that models the concentration of a neurotransmitter in the extracellular space of a brain region. This analysis yielded the first two terms in an asymptotic expansion of the concentration, which shows the first and second order effects of the many details and parameters in the problem. We developed a numerical solution method to verify our asymptotic results in some special cases.

Our asymptotic analysis assumed that the nerve varicosities are much smaller than the typical distance between them. Following [10], the distance between varicosities varies, but for serotonin there are approximately $2.6\times 10^{6}$ varicosities per cubic millimeter [60], which implies a distance between varicosities of approximately $l_{1}=7\,\mu\text{m}$. In Figure 1 of [61], varicosities are separated by about $l_{1}=20\,\mu\text{m}$ on average. The average radius of a varicosity in [60] is approximately $l_{0}=.3\,\mu\text{m}$. Hence, we estimate that the ratio $\varepsilon=l_{0}/l_{1}$ is indeed small,

[TABLE]

Of course, our model ignores various features of the actual biological system. For one, the rate of neurotransmitter release is likely not constant when a neuron is firing, and the varicosities do not perfectly absorb neurotransmitter when the neuron is not firing. In fact, a more detailed model might allow the neurotransmitter release rate and/or firing rate to depend on the extracellular concentration, and these may be interesting avenues for further research. Furthermore, we have ignored the presence of obstacles in the extracellular space which restrict the diffusive path of neurotransmitter (such as the actual neurons in the extracellular space; we have merely included the bulbous varicosities on a neuron). Indeed, the extracellular space comprises about 20% of brain volume and is a very tortuous and highly connected foamlike structure formed from the interstices of cell surfaces [62]. However, using the diffusion equation to model biochemical concentrations in the extracellular space in the brain has a long history [63] and its validity is well-established for length scales of approximately $2$ to $200\,\mu\text{m}$ [64, 65, 62]. The tortuosity of the extracellular space is often modeled by a reduction in the effective diffusion coefficient compared to diffusion in free solution [65]. The tortuosity could also potentially be modeled by a space-dependent diffusion coefficient.

Perhaps more importantly, our model does not include non-neuronal cells which can take up neurotransmitter. While the brain contains $10^{11}$ neurons, it actually contains at least ten times as many glial cells [66]. Furthermore, some glial cells are known to take up neurotransmitter [67, 68], though their neurotransmitter absorption rate is likely not as fast as that of nerve varicosities. In fact, a recent study combining experiments and mathematical modeling suggests that this uptake pathway is important for serotonin [69]. An important extension of our model and analysis would be to understand how these glia affect the extracellular neurotransmitter concentration. We expect that the steady-state neurotransmitter concentration would no longer be spatially constant if we included this uptake by glial cells. However, it remains to determine the magnitude of the spatial variation that would be introduced by this uptake.

An additional avenue for future work is investigating the variance of the neurotransmitter concentration in the extracellular space. Previous analysis of a simplified model in one space dimension showed that this variance is approximately constant outside a small neighborhood of a varicosity, but quantitative estimates are lacking for models in three space dimensions. This analysis would require significant extensions of the asymptotic and numerical methods employed in the present work, since calculating the variance of a three-dimensional stochastic model requires solving a certain boundary value problem in six space dimensions which couples at the boundaries to the mean of the stochastic model [21].

Acknowledgements. SDL thanks Janet Best, Fred Nijhout, and Michael Reed for stimulating discussions.