Synaptic Channel Modeling for DMC: Neurotransmitter Uptake and Spillover in the Tripartite Synapse
11 Synaptic Channel Modeling for DMC:Neurotransmitter Uptake and Spillover in theTripartite Synapse
Sebastian Lotter, Arman Ahmadzadeh, and Robert Schober
Abstract
In Diffusive Molecular Communication (DMC), information is transmitted by diffusing molecules. Synap-tic signaling, as a natural implementation of this paradigm, encompasses functional components that, onceunderstood, can facilitate the development of synthetic DMC systems. To unleash this potential, however, athorough understanding of the synaptic communication channel based on biophysical principles is needed.Since synaptic transmission critically depends also on non-neural cells, such understanding requires theconsideration of the so-called tripartite synapse. In this paper, we develop a comprehensive channel modelof the tripartite synapse encompassing a three-dimensional, finite-size spatial model of the synaptic cleft,molecule uptake at the presynaptic neuron and at glial cells, reversible binding to individual receptors at thepostsynaptic neuron, and spillover to the extrasynaptic space. Based on this model, we derive analytical timedomain expressions for the channel impulse response (CIR) of the synaptic DMC system and for the number ofmolecules taken up at the presynaptic neuron and at glial cells, respectively. These expressions provide insightinto the impact of macroscopic physical channel parameters on the decay rate of the CIR and the reuptakerate, and reveal fundamental limits for synaptic signal transmission induced by chemical reaction kinetics andthe channel geometry. Adapted to realistic parameters, our model produces plausible results when compared toprevious experimental and simulation studies and we provide results from particle-based computer simulationsto further validate the analytical model. The proposed comprehensive channel model admits a wide range ofsynaptic configurations making it suitable for the investigation of many practically relevant questions, such asthe impact of glial cell uptake and spillover on signal transmission in the tripartite synapse.
I. I
NTRODUCTION
In nature, information exchange between biological entities, such as cells or organs, is often based onthe release, propagation, and sensing of molecules. This process is called Molecular Communication(MC). Although traditionally studied by biologists and medical scientists, it has recently also attractedinterest from the communications research community [2]. MC is envisioned to open several exciting
This paper was accepted for presentation in part at the International Conference on Communications (ICC), 2020. a r X i v : . [ q - b i o . S C ] M a y new application areas for which communication at nano-scale is key, but can not be implementedusing traditional wireless communication systems [3]. One of these applications is the deployment andoperation of synthetic cells in the human body for the purpose of tumor detection and treatment [4];another one is the development of brain-machine interfaces (BMIs) for the detection and replacementof dysfunctional neural units [5]. In both cases, communication at cell level, either between syntheticcells or between synthetic and natural cells, is required.MC systems in which molecules propagate solely via Brownian motion, referred to as DiffusiveMolecular Communication (DMC) systems, constitute a promising candidate for synthetic MC im-plementations, as they require neither special infrastructure, nor external energy supply. The designof synthetic DMC systems, however, poses several challenges.Firstly, in environments for which synthetic DMC is intended, such as the human body, the amountof molecules available for transmission is typically limited and hence the employed communicationscheme needs to be extremely energy efficient. The concept of Energy Harvesting (EH) enablesultra-low-power applications in traditional communications [6] and several methods to leverage italso for DMC have been proposed recently [7]–[9]. The models considered in [7], [8], however, arebased on free-space propagation and can thus not exploit the channel characteristics, while in [9] anon-biological EH mechanism is considered which may not be applicable in all DMC environments.Secondly, because Brownian motion is an undirected propagation mechanism, DMC channels aretypically dispersive and thus susceptible to inter-symbol interference (ISI). This problem is well-known and, accordingly, several approaches for ISI mitigation in DMC have been proposed in theliterature, including ISI-aware modulation schemes [10], forward error-correction codes [11], [12],and equalization techniques [13]. While [10]–[13] focus mainly on the transmitter and the receiver,in [14], [15] enzymatic degradation of information molecules in the channel is considered for themitigation of ISI. This approach is interesting, because it leverages capabilities inherent to the MCchannel and does not increase the complexity of the transmitter or receiver. It does, however, incurhigher energy cost for the production of enzymes and additional signaling molecules, and is thus notnecessarily energy efficient.While the approach in [14] was inspired by the neuromuscular junction, in this paper, we considera different natural DMC system, namely molecular synaptic transmission between two neurons.Abstracted in communication terms, here, the presynaptic neuron (transmitter) encodes a sequenceof electrical impulses (data stream) into a spatio-temporal neurotransmitter release pattern (molec-ular signal) which propagates through the synaptic cleft (channel) and is finally received by the postsynaptic neuron (receiver) where the neurotransmitters activate membrane receptors. In addition,trans-membrane transporter proteins at the presynaptic neuron and at surrounding glial cells clearsignaling molecules from the synaptic space and return them to the presynaptic axon terminal forrecycling and future reuse [16]. In this way, the synaptic transmission is effectively terminated and,in addition, neurotransmitters are prevented from escaping from the synapse and causing interferenceat neighboring synapses (spillover) [17]. Depending on the specific type of synapse, the reuptake ateither the presynaptic neuron [18] or at surrounding neuroglia [16] can be dominant.It is known that the reuptake of neurotransmitters is critical for synaptic communication; several se-vere neurological diseases, including Amyotrophic Lateral Sclerosis (ALS), Alzheimer Disease (AD),Ischemia, and Epilepsy are associated with overactivation of excitatory receptors due to dysfunctionalreuptake [19], [20]. Furthermore, reuptake is also critical for the supply with neurotransmitters atthe presynaptic neuron because neurons themselves cannot synthesize neurotransmitters [21]. Hence,through the lens of a communications engineer, ISI mitigation and EH are vital for natural synaptictransmission, and thus, the study of neurotransmitter reuptake is important for synthetic MC systemdesign.From a design perspective, in turn, it is most relevant to understand how the synaptic channelparameters impact reuptake and under which conditions reliable and energy-efficient communicationis possible. To this end, a useful mathematical channel model should be physical to the extent that itcaptures the distinctive geometric features of the synaptic cleft and the nature of the most importantchemical reactions, while it should retain analytical tractability to the extent that insight can still beprovided.Most existing mathematical models for synaptic communication in the absence of glial cells,however, are based on free-space propagation [22], [23] and thus cannot capture the specifics ofthe synaptic channel. In [24], the synaptic cleft is modeled as infinite region bounded by two parallelplanes which represent the presynaptic and postsynaptic membrane, respectively. While the model in[24] incorporates presynaptic reuptake, it does not account for reversible postsynaptic binding and isalso not suited for three-dimensional analysis involving uptake by glial cells. Among the mathematicalmodels proposed for signal transmission in the tripartite synapse, see e.g. [25]–[27], to the best ofthe authors’ knowledge, none features spatial resolution.The main contribution of this paper is a comprehensive, three-dimensional channel model of thetripartite synapse encompassing a finite-size spatial model of the synaptic cleft, molecule uptakeat the presynaptic neuron and at glial cells, and reversible binding to individual receptors at the postsynaptic neuron. Furthermore, the model can be used to model spillover to the extrasynapticspace or communication in a synaptic channel without glial cells. Based on this model, we firstderive analytical time domain expressions for the channel impulse response (CIR) of the synapticDMC system and for the number of molecules taken up at the presynaptic neuron and at glial cells,respectively. Then, we use these expressions to characterize the fundamental limits that the channelgeometry and the chemical reaction kinetics impose on the decay rate of the CIR and the reuptakerate, respectively. Our analysis shows that there exists an upper bound on the CIR decay resultingfrom the channel geometry which is independent of the reaction kinetics and, similarly, it shows theexistence of an upper bound induced by the reaction kinetics which is independent of the channelgeometry. Furthermore, we also examine the impact of the individual macroscopic physical channelparameters, such as the presynaptic reuptake rate, on the asymptotic signal decay. Thus, our modelprovides a framework for understanding the fundamental trade-offs between different synaptic channelparameters limiting signal transmission in the tripartite synapse in the presence of glial cell moleculeuptake and spillover.To ensure the biological significance of our model, realistic values for the various model parametersare adopted and the model outcome is compared with results from experimental studies. Our analyticalresults are further validated by particle-based computer simulations and experimentally shown toprovide an unbiased estimator for the CIR. When applied to different glial cell configurations, ourmodel provides interesting novel insights into the potential of glial cell molecule uptake for ISImitigation.The mathematical model presented in this paper is a generalization of the model previously reportedin [1] from one to three spatial dimensions. Since the one-dimensional model in [1] does not addressthe impact of molecule uptake by glial cells and spillover, respectively, the three-dimensional modelpresented in this paper is far more general and hence allows for richer analyses.The remainder of this paper is organized as follows: In Section II, we state the system model andthe main assumptions used in Section III to derive the time- and space-dependent synaptic moleculeconcentration, the CIR, and the number of uptaken particles. In Section IV, different regimes for theCIR decay and the molecule reuptake rate are specified, and the asymptotic impact of the channelparameters is investigated. The particle-based simulator design is outlined in Section V, and numericalresults are presented in Section VI. Finally, the main findings are summarized in Section VII. Figure 1. Model synapse. Neurotransmitters enclosed in vesicles at the presynaptic neuron are released into the synaptic cleft, propagateby Brownian motion, and activate receptors at the postsynaptic neuron. Binding to postsynaptic receptors is reversible. Furthermore,particles are reuptaken (and recycled) at the presynaptic neuron. The cuboid represents an abstraction of the synaptic cleft. Its orangeand green surfaces represent the membrane of the pre- and postsynaptic neuron, respectively. The blue faces and the two faces oppositeto the blue ones which are transparent in this figure, may represent the cell membrane of astroglia encapsulating the synapse, but arein principle configurable to represent any kind of reflective, partially absorbing, or fully absorbing boundary.
II. S
YSTEM M ODEL
A. Geometry and Assumptions
The synaptic cleft is abstracted in our model as a three-dimensional rectangular cuboid (seeFig. 1). This simplification appears biologically meaningful when compared to images of hippocampalsynapses obtained by electron microscopy [28, Fig. 2(H)]. Moreover, it is analytically more tractablethan the actual non-regular shaped synaptic domain, while still retaining its characteristic features. Asimilar geometric model is used in [24], with the difference that the model considered here is boundedin all dimensions, while in [24] only one dimension is bounded.Formally, we characterize the domain of the synaptic cleft in Cartesian coordinates as follows
Ω = { ( x, y, z ) | x min ≤ x ≤ x max , y min ≤ y ≤ y max , z min ≤ z ≤ z max } , (1)and the concentration of molecules in µ m − at any time t at any location within the box definedby (1) as C ( x, y, z, t ) . To ease notation, we follow the convention to suppress the dependence of C ( x, y, z, t ) on x , y , z , and t in the notation whenever possible, i.e., we simply write C instead.Then, the presynaptic and postsynaptic membranes are located at x = x min and x = x max ,respectively. The other system boundaries, y = y min , y = y max , z = z min , and z = z max , can beconfigured depending on the context the model is used in. Here, as we are interested in reuptake inthe tripartite synapse, they will mostly represent the membranes of glial cells surrounding the synapse.A different use case, also covered by this model, is the evaluation of the number of neurotransmitters that (irreversibly) leave the synapse and possibly cause interference at other synapses, i.e., spillover[17].To derive the (dimensionless) impulse response of the synaptic channel, h ( t ) , we consider instanta-neous release of N particles at time t = t at location ( x, y, z ) = ( x , y , z ) ∈ Ω . h ( t ) is then givenas the number of particles adsorbed to the receptors of the postsynaptic membrane as a function oftime. For our analysis, we make four assumptions:A1) All adsorption and desorption coefficients as well as the diffusion coefficient are constant.A2) Particles diffuse independently of each other and receptors cannot be occupied, i.e., multipleparticles may bind to a single receptor.A3) The receptors at the postsynaptic membrane are uniformly distributed over the plane x = x max .A4) Reversible adsorption to individual receptors with intrinsic binding rate κ a in µ m µ s − andintrinsic unbinding rate κ d in µ s − can be treated equivalently as reversible adsorption to ahomogeneous surface with effective binding rate κ a in µ m µ s − and the same unbinding rate, κ d .A1) is a reasonable assumption if the time frame under consideration is sufficiently short. A2) can bejustified for low molecule concentrations [29], [30] and small-sized receptors [31]. Together with A1)it renders the system under consideration linear and time-invariant, such that the derivation of the CIRis actually meaningful. A3) is reasonable if we assume that the so-called postsynaptic density [32], thepart of the postsynaptic membrane which contains most receptors, extends over the entire postsynapticsurface under consideration. Finally, A4) is not obvious. For irreversible adsorption to an otherwisereflective surface covered by partially adsorbing disks, boundary homogenization, i.e., the assumedequivalence, has been justified in [33] and quantitatively refined with results from computer simulationin [34]. However, it is not clear if similar techniques can be applied for reversible reactions, too. Thereare two main issues to consider here. First, the steady-state fluxes are substantially different; the netflux at a reversibly adsorbing boundary is zero at steady-state, while for irreversible adsorption it isnonzero. Second, desorption alters the spatial concentration profile of particles near the boundary;particles are more concentrated near receptors compared to irreversible adsorption [35]. The first issuecan be resolved by calibrating the effective adsorption coefficient to the homogenized surface suchthat it correctly reproduces the steady-state flux to the patchy surface (see Section V), while thesecond issue can be resolved in a biologically plausible manner by assuming that particles may notre-adsorb immediately after unbinding [36].According to Fick’s second law of diffusion, (average) Brownian particle motion can be described by the following partial differential equation: ∂C∂t = D ∆ C, (2)where D > denotes the particle diffusion coefficient in µ m µ s − , and ∆ is the Laplace operator inCartesian coordinates, defined as ∆ f = ∂ f∂x + ∂ f∂y + ∂ f∂z . B. Pre- and Postsynaptic Neurons
Particle reuptake at the presynaptic neuron is modeled as irreversible adsorption to the homogeneousboundary at x = x min (possibly after appropriate boundary homogenization) with reuptake coefficient κ r in µ m µ s − . To simplify notation, we set x min = 0 and x max = l x . Then, the boundary conditionmodeling presynaptic reuptake is given as [37] D ∂C∂x = κ r C, at x = 0 , (3)where κ r ≥ . A boundary condition of the form in (3) is referred to as boundary condition of thethird kind [38].To model reversible adsorption to the postsynaptic boundary, we adopt the backreaction boundarycondition [39] and obtain D ∂C∂x = − κ a C − κ d (cid:90) t D ∂C∂x d τ, at x = l x , (4)where κ a > and κ d ≥ . Although κ a = 0 would technically be also valid, the signal at thepostsynapse would in this case be identical to , so we restrict our considerations to κ a > . C. Molecule Uptake at Glial Cell
For the boundaries at y = y min , y = y max , z = z min , and z = z max , we use boundary conditions ofthe third kind, because they provide most flexibility to our model. In particular, perfectly absorbingand perfectly reflecting surfaces are covered as special cases. For neurotransmitter uptake by glialcells, the situation is equivalent to the presynaptic membrane and thus this choice is justified by thesame argument as in Section II-B.Setting y min = 0 , y max = l y , z min = 0 , and z max = l z , this yields k ,j ∂C∂j = κ ,j C at j = 0 , j ∈ { y, z } , (5) k ,j ∂C∂j = − κ ,j C at j = l j , j ∈ { y, z } , (6)where k i,j ∈ { , D } , κ i,j ≥ in µ m µ s − is the adsorption coefficient (possibly after boundaryhomogenization) at boundary j = 0 for i = 1 and j = l j for i = 2 , and k i,j and κ i,j are not both tovanish at the same time, i.e., k i,j + κ i,j > . D. Particle Release
We consider the instantaneous release of N particles at t . Without loss of generality (w.l.o.g.), weset t = 0 and N = 1 . The former assumption is not restrictive because the system is time-invariant.The latter one is not restrictive because the system is linear and, hence, results can simply be rescaledby N .The instantaneous release of particle is then modeled as the initial value C = C ( x, y, z,
0) = δ ( x − x ) δ ( y − y ) δ ( z − z ) , ( x , y , z ) ∈ Ω , (7)where δ ( · ) denotes the Dirac delta function.In reality, when we consider multiple vesicular releases of neurotransmitters into the synaptic cleft,these do not always happen at the same release position, neither is the release instantaneous. However,as we will derive the Green’s function of the system, our result can be easily adapted to more realisticrelease models by using the theory of Green’s functions [40].
E. Channel Impulse Response and Number of Reuptaken Particles
The number of particles adsorbed to the postsynaptic boundary at time t equals the sum of allpositive and negative fluxes over this boundary in the interval [0; t ] . Thus, once the solution to (2)–(7)is found, the CIR, h ( t ) , can be obtained as the accumulated net flux over the entire postsynapticboundary, h ( t ) = t (cid:90) l y (cid:90) l z (cid:90) − D ∂C ( x, y, z, τ ) ∂x (cid:12)(cid:12)(cid:12)(cid:12) x = l x d z d y d τ. (8)For our further analysis, we denote h ( t ) in the special case of reflective boundaries in y and z , i.e., κ i,j = 0 , ∀ i ∈ { , } , j ∈ { y, z } , by h r ( t ) .Similarly, the number of reuptaken particles, i ( t ) , is obtained as the cumulative flux over thepresynaptic boundary, i ( t ) = t (cid:90) l y (cid:90) l z (cid:90) D ∂C ( x, y, z, τ ) ∂x (cid:12)(cid:12)(cid:12)(cid:12) x =0 d z d y d τ. (9)The number of reuptaken particles is a relevant quantity, because, first, it provides an upper boundon the ISI caused by the signaling molecules remaining in the synaptic cleft after transmission and,second, since neurons cannot synthesize neurotransmitters on their own [21], it provides a lowerbound on the time that the presynapse needs to recover from signal transmission, such that it canrefill the pool of its readily releasable vesicles. If the boundaries in y and z are reflective ( κ i,j = 0 , ∀ i ∈ { , } , j ∈ { y, z } ) and adsorption to thepostsynaptic boundary is actually reversible ( κ d > ), eventually all particles are reuptaken. In thisspecial case, we denote i ( t ) as i r ( t ) and have lim t →∞ i r ( t ) = 1 . (10)As D ∂C∂x (cid:12)(cid:12)(cid:12)(cid:12) x =0 = κ r C ≥ , (11)we can then interpret i r ( t ) as cumulative distribution function for the probability that a particle releasedat t = 0 is reuptaken in the interval [0; t ] .Furthermore, for the investigation of glial cell molecule uptake, we define the number of moleculesuptaken at the glial cell as i − ( t ) = t (cid:90) l x (cid:90) l z (cid:90) k ,y ∂C ( x, y, z, τ ) ∂y (cid:12)(cid:12)(cid:12)(cid:12) y =0 − k ,y ∂C ( x, y, z, τ ) ∂y (cid:12)(cid:12)(cid:12)(cid:12) y = l y d z d x d τ + t (cid:90) l x (cid:90) l y (cid:90) k ,z ∂C ( x, y, z, τ ) ∂z (cid:12)(cid:12)(cid:12)(cid:12) z =0 − k ,z ∂C ( x, y, z, τ ) ∂z (cid:12)(cid:12)(cid:12)(cid:12) z = l z d y d x d τ, (12)and the total number of uptaken particles as i tot ( t ) = i ( t ) + i − ( t ) . (13)If κ d > , as particles cannot leave the system over any other boundary, we have lim t →∞ i tot ( t ) = 1 (14)and i tot ( t ) is the cumulative distribution function for the probability that a particle is uptaken in theinterval [0; t ] at either the glial cell or the presynaptic neuron.III. A NALYTICAL C HANNEL M ODEL
A. Molecule Concentration
If our system model would comprise only boundaries of the third kind, the solution of the three-dimensional problem would be just a multiplication of solutions to the corresponding one-dimensionalproblems and could easily be obtained from the literature [40]. However, the particular challenge insolving (2)–(7) is the temporal coupling induced by the backreaction boundary condition (4). Proposition 1:
Let κ r , κ d ≥ , κ a > , k i,j ∈ { , D } , κ i,j ≥ , and k i,j + κ i,j > , ∀ i ∈ { , } , j ∈{ y, z } . Then, the unique solution to (2)–(7) is C ( x, y, z, t ) = ∞ (cid:88) m =1 ∞ (cid:88) p =1 (cid:20) κ d κ a + l x κ d l y l z { } ( κ r ) { } ( β m ) { } ( γ p )+ ∞ (cid:88) n =1 A m,pn B m C p C ( x , y , z ) C m,pn,x ( x ) C m,y ( y ) C p,z ( z ) exp( − λ m,pn t ) (cid:35) , (15)where C m,pn,x ( x ) = κ r sin( α m,pn x ) + Dα m,pn cos( α m,pn x ) , (16) C m,y ( y ) = κ ,y sin( β m y ) + k ,y β m cos( β m y ) , (17) C p,z ( z ) = κ ,z sin( γ p z ) + k ,z γ p cos( γ p z ) , (18) C ( x , y , z ) = C m,pn,x ( x ) C m,y ( y ) C p,z ( z ) , (19) A m,pn = 2 (cid:2) ( α m,pn ) D ( λ m,pn − κ d ) + ( λ m,pn ) κ a (cid:3) / { ( λ m,pn − κ d ) (cid:2) D ( α m,pn ) ( κ r ( λ m,pn − κ d ) + κ a λ m,pn ) + κ r κ a Dλ m,pn (cid:3) + 2 κ r κ a κ d D ( α m,pn ) + κ r κ a D ( λ m,pn ) + 2 κ a κ d D ( α m,pn ) + l x (cid:2) κ r D ( α m,pn ) ( λ m,pn − κ d ) + κ a D ( α m,pn ) ( λ m,pn ) + D ( α m,pn ) ( λ m,pn − κ d ) + κ r κ a ( λ m,pn ) (cid:3) } , (20) B m = (2 − { } ( β m ))( β m k ,y + κ ,y ) / (cid:2) l y ( κ ,y + k ,y β m )( κ ,y + k ,y β m )+ κ ,y k ,y ( κ ,y + k ,y β m ) + κ ,y k ,y ( κ ,y + k ,y β m ) (cid:3) , (21) C p = (2 − { } ( γ p ))( γ p k ,z + κ ,z ) / (cid:2) l z ( κ ,z + k ,z γ p )( κ ,z + k ,z γ p )+ κ ,z k ,z ( κ ,z + k ,z γ p ) + κ ,z k ,z ( κ ,z + k ,z γ p ) (cid:3) , (22) λ m,pn = D (( α m,pn ) + β m + γ p ) , (23) α m,pn are the non-zero roots of tan( αl x ) = αD (cid:2) D ( α + β m + γ p )( κ r + κ a ) − κ r κ d (cid:3) D ( α + β m + γ p )( α D − κ r κ a ) − α κ d D , (24)with either positive real part or positive imaginary part, β m and γ p are the positive, or, if κ ,y = κ ,y = 0( κ ,z = κ ,z = 0) , the real, non-negative, roots of tan( βl y ) = β ( κ ,y k ,y + κ ,y k ,y ) β k ,y k ,y − κ ,y κ ,y , (25) and tan( γl z ) = γ ( κ ,z k ,z + κ ,z k ,z ) γ k ,z k ,z − κ ,z κ ,z , (26)respectively, and { S } ( x ) denotes the indicator function, { S } ( x ) = 1 if x ∈ S , { S } ( x ) = 0 otherwise. Proof:
Please refer to the Appendix.
Remark 1:
The aforementioned coupling introduced by (4) becomes apparent in the dependenceof the α m,pn on β m and γ p . In fact, if κ d = 0 , i.e., (4) collapses to an ordinary boundary of the thirdkind, (20) and (24) become independent of β m and γ p , and C ( x, y, z, t ) factorizes into a product ofone-dimensional solutions. Remark 2:
Although imaginary values are allowed for α m,pn , C ( x, y, z, t ) is always a real-valuedfunction, because the trigonometric functions occurring in (15) when evaluated at purely imaginaryarguments equal their hyperbolic counterparts evaluated at real arguments and the leftover imaginaryparts cancel each other out. Remark 3: C ( x, y, z, t ) is called the Green’s function of the boundary value problem (2)–(6) [40].We index the sets { β m } and { γ p } ascendingly, such that β and γ are the smallest admissiblesolutions of (25) and (26), respectively. Furthermore, we remark that (24) admits at most one imaginarysolution with positive imaginary part and index the real roots of (24) ascendingly, starting with index if (24) does not admit an imaginary root and starting with index otherwise. In the latter case, weindex the imaginary root of (24) with index . B. CIR and Number of Reuptaken Particles
Next, we use Proposition 1 to derive the CIR and the number of reuptaken particles.
1) CIR:
First, we consider the CIR.
Corollary 1:
The CIR, h ( t ) , as defined in (8) can be obtained from (15) as h ( t ) = ∞ (cid:88) n =1 ∞ (cid:88) m =1 ∞ (cid:88) p =1 A m,pn B m C p C ( x , y , z )( κ r α m,pn cos( α m,pn l x ) − D ( α m,pn ) sin( α m,pn l x )) × (cid:20) k ,y sin( β m l y ) + κ ,y β m (1 − cos( β m l y )) (cid:21) (cid:20) k ,z sin( γ p l z ) + κ ,z γ p (1 − cos( γ p l z )) (cid:21) × exp( − λ m,pn t ) − α m,pn ) + β m + γ p . (27) Proof:
Eq. (27) follows by elementary differentiation and integration, cf. (8), (15). Here, to avoid the obvious difficulties if β m = 0 or γ p = 0 in (27), we observe that this only happensif κ ,y = κ ,y = 0 or κ ,z = κ ,z = 0 and in this case lim β m → B m ( κ ,y sin( β m y ) + k ,y β m cos( β m y ))( k ,y sin( β m l y ) + κ ,y β m (1 − cos( β m l y ))) =lim β m → l y β m k ,y k ,y β m cos( β m y ) k ,y sin( β m l y ) = 1 (28)and, similarly, lim γ p → C p ( κ ,z sin( γ p z ) + k ,z γ p cos( γ p z ))( k ,z sin( γ p l z ) + κ ,z γ p (1 − cos( γ p l z ))) =lim γ p → l z γ p k ,z k ,z γ p cos( γ p z ) k ,z sin( γ p l z ) = 1 (29)and thus, if β m = 0 or γ p = 0 , we simply replace the corresponding expressions in (27) with theselimits. This argument could be made more rigorous by considering the constant terms correspondingto β m = 0 , γ p = 0 separately in (27), but, here, we favor a unified treatment over mathematical rigor. Corollary 2:
If the boundaries in y and z are reflective, i.e., κ ,y = κ ,y = κ ,z = κ ,z = 0 , (27)simplifies to h r ( t ) = ∞ (cid:88) n =1 A , n C , n,x ( x ) (cid:0) κ r cos( α , n l x ) − Dα , n sin( α , n l x ) (cid:1) exp( − D ( α , n ) t ) − α , n , (30)where { β m } and { γ p } are indexed such that β = γ = 0 .If κ r > , κ d > , (30) simplifies further to h r ( t ) = ∞ (cid:88) n =1 A , n C , n,x ( x )( κ r cos( α , n l x ) − Dα , n sin( α , n l x )) exp( − D ( α , n ) t ) α , n . (31) Proof: If κ ,y = κ ,y = κ ,z = κ ,z = 0 , (25) and (26) simplify to sin( βl y ) = 0 and sin( γl z ) = 0 ,respectively, and all terms corresponding to non-zero roots β m , γ p in (27) vanish. Then, using thesubstitutions from (28) and (29), (30) follows directly from (27).Eq. (31) follows from the observation that, if κ r > , κ d > , all released particles are eventuallyreuptaken and, hence, lim t →∞ h r ( t ) = 0 .Eq. (30) is exactly the result that was obtained for the one-dimensional model in [1].Since (31) is a weighted sum of decaying exponentials, we expect that h r is dominated by theexponential with the slowest decay for large values of t . Thus, we obtain the following approximationfor large t : h r ( t ) ≈ A , C , ,x ( x ) (cid:0) κ r cos( α , l x ) − Dα , sin( α , l x ) (cid:1) exp( − D (cid:0) α , (cid:1) t ) α , (32)The accuracy of this approximation will be investigated in Section VI, cf. Fig. 3. For non-reflective boundaries, h ( t ) is asymptotically still dominated by the first summand in (27),and the first-term approximation h ( t ) ≈ A , B C C ( x , y , z ) (cid:104) κ r α , cos( α , l x ) − D (cid:0) α , (cid:1) sin( α , l x ) (cid:105) × (cid:20) k ,y sin( β l y ) + κ ,y β (1 − cos( β l y )) (cid:21) (cid:20) k ,z sin( γ l z ) + κ ,z γ (1 − cos( γ l z )) (cid:21) × exp( − λ , t ) (cid:0) α , (cid:1) + β + γ (33)provides an accurate estimate of the CIR decay, cf. Fig. 5.This observation concludes the derivation of the CIR.
2) Number of Reuptaken Particles:
Now, we consider the number of particles taken up at thepresynaptic neuron and at the glial cell.
Corollary 3:
The number of reuptaken particles at the presynaptic neuron as defined in (9) is givenas follows i ( t ) = κ r ∞ (cid:88) n =1 ∞ (cid:88) m =1 ∞ (cid:88) p =1 α m,pn A m,pn B m C p C ( x , y , z ) (cid:20) k ,y sin( β m l y ) + κ ,y β m (1 − cos( β m l y )) (cid:21) × (cid:20) k ,z sin( γ p l z ) + κ ,z γ p (1 − cos( γ p l z )) (cid:21) − exp( − λ m,pn t )( α m,pn ) + β m + γ p . (34)The number of particles uptaken at the glial cell, i − ( t ) as defined in (12) is given by i − ( t ) = ∞ (cid:88) n =1 ∞ (cid:88) m =1 ∞ (cid:88) p =1 A m,pn B m C p C ( x , y , z ) (cid:20) D sin( α m,pn l x ) + κ r α m,pn (1 − cos( α m,pn l x )) (cid:21) × (cid:26) [ κ ,y C m,y (0) + κ ,y C m,y ( l y )] (cid:20) k ,z sin( γ p l z ) + κ ,z γ p (1 − cos( γ p l z )) (cid:21) + (cid:20) k ,y sin( β m l y ) + κ ,y β m (1 − cos( β m l y )) (cid:21) [ κ ,z C p,z (0) + κ ,z C p,z ( l z )] (cid:27) × − exp( − λ m,pn t ) λ m,pn . (35) Proof:
We plug (3) into (9) and apply elementary integration to (15) to obtain (34). Similarly,we use (5) and (6) to obtain (35).
Corollary 4:
If the boundaries in y and z are reflective, i.e., κ ,y = κ ,y = κ ,z = κ ,z = 0 , (34)simplifies to i r ( t ) = κ r ∞ (cid:88) n =1 A , n C , n,x ( x ) 1 − exp( − D ( α , n ) t ) α , n , (36)where, again, we have set β = γ = 0 . If κ d > , (36) simplifies to i r ( t ) = 1 − κ r ∞ (cid:88) n =1 A , n C , n,x ( x ) exp( − D ( α , n ) t ) α , n . (37) Proof:
We use the substitutions from (28) and (29) to obtain (36) from (34). If κ d > , we canfurther exploit (10) to obtain t →∞ i r ( t ) = κ r ∞ (cid:88) n =1 A , n C , n,x ( x ) 1 α , n , (38)and (37) follows.Hence, we note that, if the boundaries in y and z are reflective, the three-dimensional domain Ω is in terms of h ( t ) and i ( t ) equivalent to the one-dimensional domain [0; l x ] .From (37), we obtain i r ( t ) ≈ − κ r C , ,x ( x ) A , exp( − D (cid:0) α , (cid:1) t ) α , =: i r ( t ) . (39)as an approximation for the number of particles taken up at the presynaptic neuron. The accuracy ofthis approximation is investigated in Section VI, cf. Fig. 2.IV. CIR D ECAY R ATE
Now, in order to develop a better understanding of the limits and mechanisms that determine thepostsynaptic signal as well as the reuptake of molecules, we first characterize α , for the special caseof reflective boundaries in y and z , and then discuss the more general case of partially absorbingboundaries. A. Reflective Boundaries
From (32) and (39), we observe that, if the boundaries in y and z are reflective, the smallest root of(24), α , , dictates the decay rate of the CIR, and, at the same time, governs the rate of the presynapticreuptake. For ease of notation, we drop the dependence of α , on β and γ and write α instead.Now, to characterize α , we consider both sides of (24), which for κ ,y = κ ,y = κ ,z = κ ,z = 0 simplify to tan( αl x ) = α D ( κ r + κ a ) − κ r κ d α [ α D − ( κ r κ a + κ d D )] =: u ( α ) . (40)All solutions of this equation are real and we seek the smallest positive one. Proposition 2:
Let t l x = l x D denote the time that a particle needs (on average) to travel distance l x if it diffuses with diffusion coefficient D . Then, the following inequality holds for the rate whichdetermines the exponential decay in (32) and (39), Dα : Dα < min (cid:26) π t l x , κ d κ r κ r + κ a (cid:27) (41) Proof:
Let us first examine the middle term in (40), which we denote by u ( α ) . It has two non-negative poles, one at α = 0 and one at α = √ κ r κ a + κ d D/D =: α (2) . It can easily be shown that lim α → + u ( α ) = + ∞ and lim α → α − (2) u ( α ) = −∞ . Furthermore, u ( α ) possesses exactly one positiveroot at α = κ r κ d D ( κ r + κ a ) . The left-hand side of (40) is periodic with periodicity πl x and has poles at π l x + kπl x and roots at kπl x , k ∈ Z . Within each interval ( kπl x − π l x ; kπl x + π l x ) , tan( αl x ) is strictly monotonicallyincreasing and continuous and thus takes on all values from −∞ to + ∞ . Hence, we conclude thatthe smallest root of (40), α , (i) is contained in the interval (0; π l x ) , (ii) is smaller than (cid:113) κ r κ d D ( κ r + κ a ) .Now, from (i), we obtain Dα < Dπ l x = π t l x , (42)and (ii), on the other hand, yields Dα < κ d κ r κ r + κ a . (43)This completes the proof.Thus, because κ r κ r + κ a < the decay of h r ( t ) can never be faster than the rate of the desorptionreaction at the postsynapse, κ d . On the other hand, channel clearance is fast, if presynaptic reuptakedominates postsynaptic binding, i.e., κ r (cid:29) κ a , and slow, if molecules tend to bind to the postsynapseinstead of being reuptaken, i.e., κ a (cid:29) κ r .From (42) and (43), we see that the decay of the CIR can either be limited by the time thatparticles need to cross the synapse or by the reaction kinetics at the pre- and postsynapse, dependingon whether π t lx < κ d κ r κ r + κ a or κ d κ r κ r + κ a < π t lx . We call the regime π t lx < κ d κ r κ r + κ a diffusion-limited (in x ) , and the regime κ d κ r κ r + κ a < π t lx reaction-limited (in x ) . Thus, if the two neuronsare very close, i.e., l x is small, synaptic transmission happens in the reaction-limited regime, whilefast molecule dissociation at the postsynapse and fast presynaptic uptake favor the diffusion-limitedregime. Eqs. (42) and (43) reveal that both, diffusion and reaction, impose fundamental bounds onsignal decay and particle reuptake and, more importantly, these bounds are independent of each other,i.e., in the reaction-limited regime, independent of the geometry of the synaptic cleft, signal decayand reuptake can never exceed (43), and, similarly, in the diffusion-limited regime, even for arbitrarilyfast reactions, (42) applies and bounds the maximum achievable signal decay.So far, we have seen how the system parameters influence the (un)achievable range of α . Now,we examine their direct impact on the value of α . To this end, we fix any α ∗ ∈ (0; π l x ) and consider u instead as a function of any channel parameter which is currently under inspection, e.g. u ( κ r ) .Similarly, we consider α = α ( κ r ) as a function of that parameter. Finally, for the investigation of l x , we consider tan( αl x ) as a function of l x . Proposition 3: α is increasing in κ r and κ d and decreasing in κ a and l x , i.e.,d α d κ r , d α d κ d > , d α d κ a , d α d l x < . (44) Proof:
Due to the strict monotonicity of tan( · ) , d u ( κ r ) d κ r > implies d α d κ r > . The same holdsfor κ a and κ d , and for l x , d tan( α ∗ l x ) d l x > implies d α d l x < . By simple differentiation, we obtain d u ( κ r ) d κ r , d u ( κ d ) d κ d , d tan( α ∗ l x ) d l x > , d u ( κ a ) d κ a < and, hence, (44) follows. This completes the proof.Thus, increasing reuptake and desorption accelerates channel clearance, while larger postsynapticadsorption and synaptic cleft width slow it down. B. Non-Reflective Boundaries
In the general case of partially absorbing boundaries in y and z , the decay of the CIR is governedby the three-dimensional decay rate λ , = D (cid:104)(cid:0) α , (cid:1) + β + γ (cid:105) , see (33).Let us first consider α , . If α , is real, the discussion parallels the discussion of α with minormodifications due to β and γ . Proposition 4: If α , ∈ R , the exponential decay rate λ , of the first-term approximation of h ( t ) as defined in (33) is bounded by the following inequality: λ , < min (cid:26) π t l x + D ( β + γ ) , κ d κ r κ r + κ a (cid:27) (45) Proof:
With similar arguments as in the proof of Proposition 2, one can show that, if α , ∈ R , α , , (i) is contained in the interval (0; π l x ) , (ii) is smaller than (cid:113) κ r κ d D ( κ r + κ a ) − ( β + γ ) . Hence, weconclude that for non-reflective boundaries the upper bound on the maximum decay rate induced bythe channel geometry is λ , < π t l x + D ( β + γ ) , (46)and the reaction rate induced upper bound is the same as in the one-dimensional scenario, namely λ , < κ d κ r κ r + κ a . (47)This completes the proof.For many practical scenarios, β and γ are too small relative to π t lx so as to significantlyincrease the bound in (46). This means that for such scenarios the maximum CIR decay of a synaptictransmission system with glial cell uptake can also be achieved by an equivalent one-dimensionalsystem with high presynaptic reuptake and postsynaptic desorption rate. The exception to this ruleare synaptic transmission systems in which the channel is “wider than high”, i.e., l y and l z are relativelysmall compared to l x . Such systems can for example be found in the retina [41]. In this case, as the values of β and γ scale with the reciprocal of l y and l z , respectively, β and γ are in a range inwhich (46) is significantly larger than (42).Now, let us consider the case α , (cid:54)∈ R . In order for (24) to allow for an imaginary root, twoconditions must be fulfilled. First, there must be a flow of molecules back from the postsynapticboundary into the system, i.e., κ d > , and, second, absorption in y and/or z must be strong, i.e., β m + γ p must be large. Under these conditions, a pair of complex conjugate, purely imaginary, rootsemerges as solution of (24). Let us write α , = jθ . Then, we have λ , = D ( − θ + β + γ ) andsee that α , contributes negatively to the decay of the CIR, i.e., particle transmission in x direction slows down the signal decay.At first glance, this seems to be counter-intuitive, but actually it is a very sensible insight: First,being adsorbed to the postsynaptic membrane prevents a particle from being absorbed by the “hungry”glial cell and, second, a particle that has desorbed from the postsynaptic membrane may re-adsorbagain and in this way contribute positively to the postsynaptic signal. This is true in general, but when β m + γ p is large, uptake at glial cells is so fast that desorption from the postsynaptic membrane actseffectively as a molecule source. The signal decay in the synaptic cleft is then clearly dictated by theglial cell.Now, to complete our discussion of λ , , let us examine β , the same analysis holds for γ . Wehave already seen that if κ ,y = κ ,y = 0 , β = 0 , i.e., there is zero contribution from β to λ , .Thus, consider the case in which at least one of κ ,y and κ ,y is non-zero. Then, β is defined as thesmallest non-zero solution of (25). Proposition 5: If κ ,y + κ ,y > , the following inequalities hold: β (1) < β < π l y if β (1) < π l y (48) π l y < β ≤ πl y if β (1) > πl y (49) π l y < β < β (1) otherwise , (50)where β (1) = (cid:113) κ ,y κ ,y D and equality in (49) is attained if and only if k ,y = k ,y = 0 . Proof:
The analysis here is a bit different from the above analysis of α . First, we considera partially adsorbing boundary, k ,y , k ,y = D (cid:54) = 0 . We note that the right-hand side of (25),which we denote by v ( β ) = βD ( κ ,y + κ ,y ) β D − κ ,y κ ,y , has a root at β = 0 and a pole at β (1) = (cid:113) κ ,y κ ,y D .Furthermore, lim β → β − (1) v ( β ) = −∞ and v ( β ) < , ∀ β ∈ (0; β (1) ) . As the left-hand side of (25) isperiodic with periodicity πl y , has poles at π l y + kπl y , and is strictly monotonically increasing in eachinterval (cid:16) π l y + kπl y ; π l y + ( k +1) πl y (cid:17) , k ∈ Z , inequalities (48)–(50) follow. Now, if k ,y = k ,y = 0 , (25) simplifies to sin( βl y ) = 0 , and thus, β m = m πl y , m ∈ N . Hence, we have β = πl y and the upper boundin (49) is attained. This completes the proof.To simplify the interpretation of this result, let us consider the special case κ ,y = κ ,y = κ y .Then, β (1) > πl y if and only if κ y > π l y t ly , where t l y = l y D denotes the time that a particle needs(on average) to travel a distance l y with diffusion coefficient D . Thus, if the reaction kinetics at theboundaries are sufficiently fast relative to the particle velocity, (49) holds, and β is constrained bythe channel geometry only. Similar to before, we call this regime diffusion-limited (in y ) , but, here,the particle velocity (relative to the channel length in y ) provides not only an upper, but also a lowerbound for the contribution of particle loss in y to the overall decay rate. This implies that, if l y decreases, e.g. due to a change in the morphology of a glial cell, as long as reactions at the glialcell happen fast enough, i.e., β (1) > πl y holds, we have a stronger guarantee on the minimum signaldecay. For fixed reaction rates, however, if l y grows beyond limits, (49) holds, β approaches , andwe recover the reflective scenario.On the other hand, if the reaction kinetics are relatively slow, i.e., κ y < π l y t ly , and either of (48) or(50) applies, β is bounded by the reaction kinetics, either from below or from above. We call thisregime again reaction-limited (in y ) , though we note that in x reaction kinetics exclusively providedan upper bound .The next proposition concludes this analysis with results on the dependence of β on κ ,y , κ ,y ,and l y . Proposition 6: β as a function of κ ,y , κ ,y , and l y is increasing in κ ,y and κ ,y and decreasing in l y , i.e., d β d κ ,y , d β d κ ,y > , d β d l y < . (51) Proof:
Similar to the proof of Proposition 3, (51) follows from the monotonicity of tan( βl y ) andelementary differentiation of v ( β ) as defined in the proof of Proposition 5.V. P ARTICLE - BASED S IMULATION
To verify the analytical expressions derived in Section III, three-dimensional particle-based com-puter simulations were conducted. To this end, we adopted the simulator design from [42].
A. Simulator Design
Here, Brownian particle motion is simulated by updating the position of each point-like particleat each time step with a three-dimensional jointly independent Gaussian random vector [ X, Y, Z ] ∼ N ( × , σ I × ) , where N ( µ , C ) denotes the multivariate Gaussian distribution with mean vector µ and covariance matrix C , and × M and I M × M denote the × M all-zero vector and the M × M identity matrix, respectively. The variance σ is a function of the particle diffusion coefficient D andthe simulation time step ∆ t , σ = 2 D ∆ t . Accordingly, the root mean step length (rms) of a simulatedparticle is given by s = √ D ∆ t. (52)In the simulation, the probability that a particle is adsorbed when hitting a homogeneous boundaryis computed from the respective adsorption coefficient using [42, eq. (21)]. If a particle hits aboundary and is not absorbed, it is reflected using ballistic reflection. Particles hitting a receptor atthe postsynaptic boundary are absorbed with probabilities computed as in [42, eqs. (37), (32)] fromthe intrinsic binding rate of molecules to receptors, κ a , and the intrinsic desorption rate constant, κ d .Desorbed particles are placed apart from the boundary according to [42, eq. (35)]. B. Boundary Homogenization
In order to compare simulation data and analytical results, boundary homogenization at the post-synaptic boundary is performed, i.e., the heterogeneous surface covered by individual receptors assimulated in the particle-based simulator is identified with an equivalent homogeneous surface as mod-eled in (4). For irreversible molecule binding to individual receptors, this technique is well-developed[33], [34], [43], [44]. However, for reversible reactions, to the best of the authors’ knowledge, it isnot available yet. Indeed, the existing general theory for reversible diffusion-influenced reactions israther inaccessible for the non-expert user [45], [46]. Here, there are two main reasons why existingresults for irreversible ligand-receptor binding, such as [34], cannot be applied.First, directly after unbinding from a receptor, molecules are physically still close to this receptorand hence, the concentration of molecules at the receptor-covered surface is non-uniform. This impliesthat the probability of hitting again a receptor when hitting the surface is higher for molecules thathave just desorbed than for molecules which were never bound to a receptor. In other words, reversiblybinding molecules have a (short-term) memory, a state, which impacts their binding properties. Thisobservation is in contrast to memoryless irreversible binding. We can, however, circumvent thisproblem in a biologically plausible manner [36] by making the assumption that there is a (short-term) mechanism which prevents ligands from immediate rebinding to the same receptor. Thus, wechoose ∆ t and the receptor radius r such that s is larger than r . As the displacement of particles afterdesorption scales with s (see [42, eq. (35)]), this has the effect that desorbed particles are more likely to see a representative part of the boundary before possibly adsorbing to a receptor again. This is infact equivalent to blocking the desorbed particle for some time from rebinding and, hence, rendersmolecule binding a memoryless process.Second, boundary homogenization for irreversible reactions is performed such that the averageparticle lifetimes in presence of the receptor-covered surface and the homogeneous surface, respec-tively, are identical [34]. This requirement, however, is not meaningful for reversible reactions becausegenerally, in the absence of any irreversible molecule degradation mechanism, the lifetime of moleculesin the presence of a reversibly absorbing surface is infinite. Thus, to compute the effective adsorptioncoefficient for the homogenized boundary, κ a , from κ a , we adopt the approach from [34] and performcomputer-assisted boundary homogenization, but with the difference, that we do not require analyticaland numerical results to produce the same average particle lifetime , but instead demand matching steady-state concentrations . To this end, we conduct particle-based simulations in Ω (1) withoutpresynaptic reuptake and with reflective boundaries in y and z until the system approaches its steady-state. Next, we use our analytical result in (15) to compute the steady-state number of particlesadsorbed to the postsynaptic boundary for this scenario, h ∞ = 1 − l x κ d κ a + l x κ d = κ a κ a + l x κ d . (53)Now, we denote the simulated steady-state number of adsorbed particles as h ∞ sim and require h ∞ sim = κ a / ( κ a + l x κ d ) , which we solve for κ a to obtain κ a = l x κ d h ∞ sim / (1 − h ∞ sim ) . Next, we compute h ∞ sim formany different parameter sets and find in this way that κ a = 0 . ρκ a , (54)where ρ denotes the fraction of the postsynaptic surface covered by receptors, provides an excellentapproximation for the range of parameters that we consider (see Section VI, cf. Table I).In contrast to [34] and earlier results, (54) indicates that, for reversible binding, κ a does not varywith the receptor radius r , as long as ρ is kept fixed. This is intuitive; it is known that the receptorradius impacts the average particle lifetime (given a constant surface coverage) [47], [48]. However,this is a transient phenomenon and does not apply to the steady-state concentration considered here.VI. N UMERICAL R ESULTS
In this section, we evaluate the analytical expressions for the CIR and the number of reuptakenparticles derived in Section III and compare them with results from particle-based simulation asoutlined in Section V. For the computation of the analytical results, the infinite sums in (27), (31), (34), (35), and (37) were truncated after terms in x and terms in y and z , respectively, i.e., n ∈ { , . . . , } , m ∈ { , . . . , } , and p ∈ { , . . . , } . Furthermore, all analytical results werescaled with the number of particles released in the particle-based simulations, N . The results fromparticle-based simulation were averaged over realizations. Unless indicated otherwise, particleswere released at ( x , y , z ) = ( s, l y / , l z / , where s is defined in (52), i.e., particle release is centeredin the y - z -plane and occurs one average simulation step length from the presynaptic boundary in x .Although a release at x = 0 . would in principle also be covered by our analytical result, we chose x = s to account for the fact that also in real synapses, due to the exocytotic nature of vesicularparticle release, upon release, particles are located in front of the presynaptic membrane and not on the presynaptic membrane [49].The remainder of this section is organized as follows. In Subsection VI-A, we provide details onthe choice of the default parameters. Next, we compare in Subsection VI-B the impact of presynapticreuptake and postsynaptic desorption on the CIR and on the number of reuptaken particles, respec-tively, without uptake at the glial cell. In Subsection VI-C, molecule uptake at the glial cell subjectto different glial cell morphologies and different neurotransmitter transporter densities is considered.Finally, we conclude this section investigating the impact of the molecule release location on the CIRin the presence of spillover in Subsection VI-D. A. Choice of Default Parameters
Whenever available, sensible default values for the model parameters were adopted from experimen-tal studies of glutamatergic excitatory hippocampal synapses as indicated in Table I. Some parametersthat were not readily available were determined as follows. The binding and unbinding rates ofour simplified two-state postsynaptic reaction scheme, κ a and κ d , were computed from the bindingrates to individual postsynaptic receptors [50] and from the postsynaptic receptor density reported in[51]. The glial cell uptake rate constants κ ,y , κ ,y , κ ,z , and κ ,z were obtained from the bindingrate constant to an individual neurotransmitter transporter and the transporter density at glial cellsreported in [52] and [53], respectively. In lack of reliable experimental values for the presynapticreuptake rate constant, we assumed in accordance with [50] a three-fold larger “uptake capacity” atglial cells compared to neurons. As the “uptake capacity” is defined as the product of transport rateper individual transporter, surface to volume ratio, and transporter density at the particular membrane,assuming that the transport rate is the same at glial cells and neurons, we compute that the transporterdensity at glial cells is approximately times larger than at neural cells and in this way obtain κ r = κ ,y /
20 = 1 . × − µ m µ s − . Table IS
IMULATION PARAMETERS FOR PARTICLE - BASED SIMULATION .Parameter Default Value Reference(s)
D, k ,y , k ,y , k ,z , k ,z . × − µ m µ s − [17] N [54] l x . µ m [55] l y , l z . µ m [55] κ r . × − µ m µ s − [50] κ a . × − µ m µ s − [50], [51] κ d . × − µ s − [50], [51] κ ,y , κ ,y , κ ,z , κ ,z . × − µ m µ s − [52], [53] ∆ t . µ s r . ρ . Finally, the simulation time step, ∆ t , the radius of the model receptors, r , and the postsynapticreceptor coverage, ρ , were chosen such that A4) from Section II is fulfilled, i.e., such that boundaryhomogenization can be applied.By this choice of the model parameters, we ensure that the model outcome is qualitatively com-parable with experimentally obtained data and, indeed, despite several simplifying assumptions, thevariation of the postsynaptic receptor activation over time predicted by our model, as can be observede.g. in Fig. 5, is in very good agreement with experimental results [54] and data from previoussimulation studies [17, Fig. 3(F)]. Note that, however, we model the postsynaptic signal processingonly up to the binding (and unbinding) of molecules to receptors and do not attempt to encompassfurther processing steps that would ultimately lead to the generation of an experimentally detectableelectrical downstream signal. Thus, the comparison with experimental data cannot be done strictlyquantitatively.In the default parameter setting, α , is imaginary and β and γ are close to the reaction-inducedlower bound from (48). Hence, the system is dominated by the reaction kinetics at the glial cell. B. Presynaptic Reuptake and Postsynaptic Desorption for Reflective Case
In Figs. 2 and 3, the impact of presynaptic reuptake and postsynaptic desorption on the CIR and onthe number of presynaptically reuptaken particles, respectively, in a synapse with reflective boundariesin y and z is shown.From Fig. 2, we observe that the time required by the presynaptic neuron to reuptake a givennumber of released particles decreases with increasing reuptake rate κ r . The desorption rate κ d onthe other hand does not impact the presynaptic reuptake significantly. N u m b e r o f R e up t a k e n M o l ec u l e s i r ( t ) [ - ] κ r = 5 . e − κ r = 1 . e − κ r = 2 . e − κ d = 2 . e − κ d = 8 . e −
05 10 − − − − − − A b s o l u t e A pp r o x i m a t i o n E rr o r | i r ( t ) − i r ( t ) | [ - ] Figure 2. Left y -axis: Presynaptic reuptake in reflective case fordifferent reuptake rates κ r in µ m µ s − and different desorptionrates κ d in µ s − . The results from the particle-based simulationsare depicted as diamond markers. Right y -axis: The absolute errorof the first-term approximation (39) is shown as dashed lines. C h a nn e l I m pu l s e R e s p o n s e h r ( t ) [ - ] κ d = 2 . e − κ d = 8 . e − κ r = 1 . e − κ r = 5 . e − κ r = 2 . e − Figure 3. CIR for different presynaptic reuptake rates κ r inµ m µ s − and different desorption rates κ d in µ s − in the reflectivecase. The first-term approximation as defined in (32) is shown asdashed lines. The results from the particle-based simulations areshown as diamond markers. Fig. 3 shows how the presynaptic reuptake and the postsynaptic desorption shape the postsynapticsignal. As the presynaptic reuptake κ r increases, the CIR decays faster and its peak value decreasesslightly. This observation underlines the role of molecule reuptake for ISI mitigation. Increasing thedesorption rate κ d , on the other hand, leads to a significant decrease of the postsynaptic signal strength.Hence, we conclude that presynaptic reuptake impacts both, the presynaptic and postsynaptic sides,while the effect of desorption is only local, at the postsynaptic side.We further observe in Figs. 2 and 3 that the first-term approximation as proposed in (39) and (33),respectively, provides an excellent approximation to quickly evaluate the reuptake dynamics and theCIR decay. For the presynaptic reuptake, the approximation is accurate even for very small valuesof t . In fact, we observed a significant decrease in the approximation quality only when pathologicalscenarios, such as a very small cleft width, were considered.For the parameters considered here, κ d κ r κ r + κ a (cid:28) π t lx and, according to Subsection IV-A and (41),the CIR decay is limited by the reaction kinetics. Note that, although in Fig. 3 the signal seems todecay slower when κ d is increased, actually only the constant in front of the exponential in (32)decreases, but the decay rate of the exponential itself increases in accordance with the results fromSection IV. C. Molecule Uptake at Glial Cells
In this section, we explore how the morphology and the reaction kinetics of the glial cell impactthe CIR and the total number of uptaken particles. This question is of practical importance because both parameters are variable in nature and have been shown to impact the synaptic channel clearance[56].We consider the default set of parameters as listed in Table I as baseline and then increase thedistance to the glial cell, reflected in our model by l y and l z , first by a factor of . to l y = l z = 0 . µ m and then by a factor of to l y = l z = 0 . µ m . Similarly, we increase the glial cell uptake rate constantfirst by a factor of . to κ ,y = κ ,y = κ ,z = κ ,z = 3 . × − µ m µ s − and then by a factor of to κ ,y = κ ,y = κ ,z = κ ,z = 5 . × − µ m µ s − .The results are shown in Figs. 4 and 5. We observe that the effects of increasing the distance fromthe transmitter release site to the glial cell and increasing the uptake rate at the glial cell on the CIRare opposite. First, let us compare curves corresponding to low, medium and large distance to theglial cell in Fig. 5. Clearly, the further the distance from release site to glial cell, the larger is thepeak value of the CIR and the slower it decays. This is intuitive, because, as the distance to the glialcell membrane increases, more signaling molecules can reach the postsynaptic membrane before theyare taken up by the glial cell. On the other hand, comparing the results for low, medium, and highuptake rates in Fig. 5, we observe that increasing the uptake rate clearly shortens the CIR and, atthe same time results in a smaller peak value. Fig. 4 confirms that the different CIR shapes actuallyhave their correspondence in the shapes of the molecule uptake. Here, larger distances lead to slower uptake, while higher uptake rates lead to faster uptake.Now, we consider the effect of the simultaneous change of parameters. First, we observe from Fig. 5that the simultaneous increase of distance and uptake rate by the same factor does not impact the CIRshape, i.e., the curves corresponding to ( l y,z , κ ) = { (0 . , . × − ) , (0 . , . × − ) , (0 . , . × − ) } , where κ = κ ,y = κ ,y = κ ,z = κ ,z , practically coincide. The same applies to the numberof uptaken particles as shown in Fig. 4. Thus, in the regime we are considering, the distance to theglial cell and the glial cell uptake rate can be traded-off one for the other.On the other hand, if we compare the curves in Fig. 5 corresponding to ( l y,z , κ ) = { (0 . , . × − ) , (0 . , . × − ) } , we observe that increasing the glial cell uptake rate can also overcompensate for an increased distance to the glial cell. Here, the distance is increased by a factor of . , while theuptake rate is increased by a factor of . relative to the default values. Most interestingly, we observethat the peak value of the CIR remains unchanged while the CIR decays faster. This is because whenincreasing the distance, molecules will typically hit the glial cell membrane later , but will then betaken up faster . This finding suggests that strong molecule uptake by the glial cell provides a potentialmechanism for selectively suppressing the undesired parts of the postsynaptic signal and hence for N u m b e r o f T o t a l R e up t a k e n M o l ec u l e s i t o t ( t ) [ - ] l { y,z } = 0 . l { y,z } = 0 . l { y,z } = 0 . κ = 2 . e − κ = 3 . e − κ = 5 . e − Figure 4. Total number of reuptaken particles for differentdistances to the glial cell l y , l z in µ m , and different glial cell uptakerates κ ,y = κ ,y = κ ,z = κ ,z = κ in µ m µ s − . The resultsfrom the particle-based simulations are depicted by the diamondmarkers. C h a nn e l I m pu l s e R e s p o n s e h ( t ) [ - ] l { y,z } = 0 . l { y,z } = 0 . l { y,z } = 0 . κ = 2 . e − κ = 3 . e − κ = 5 . e − Figure 5. CIR for different distances to the glial cell l y , l z in µ m ,and different glial cell uptake rates κ ,y = κ ,y = κ ,z = κ ,z = κ in µ m µ s − . The results from the particle-based simulations aredepicted by the diamond markers. The first-term approximation(33) is shown as dashed lines. ISI mitigation.Fig. 5 shows also that the first-term approximation as defined in (33) provides an accurate ap-proximation to h ( t ) shortly after the CIR has peaked. For t > . it is not distinguishable from h ( t ) .A natural question arising from the preceding discussion is, to which extent increasing the glial celluptake rate can speed up the CIR decay. But here, the bound for β derived in Proposition 5 comesin; increasing the glial cell uptake rate has only an effect on the CIR decay as long as the diffusion-induced upper bound from (49) is not reached. Hence, if the glial cells are far apart, increasing theuptake rate does not significantly increase the CIR decay. D. Vesicle Release Location and Spillover
As explained in [19, Section 3.2.6], the ability of glial cells to mitigate intersynaptic cross-talk,i.e., co-channel interference, which is due to the escape of signaling molecules from the synapse andtheir diffusive migration to neighboring synapses, highly depends on the position of the glial celltransporters relative to the molecule release site. Here, we investigate this phenomenon by modelingspillover to neighboring synapses as particle flux over the boundaries z = 0 and z = l z . In particular,we make the assumption that particles do not come back to the synapse once they have escaped fromit and model spillover in this way as irreversible perfect adsorption. Hence, we set k ,z = k ,z = 0 tomake the boundaries in z perfectly adsorbing. In reality, the distance from a particular synapse to neighboring synapses depends highly on thetype of synapse under consideration and even within the same type is subject to variations by morethan a factor of [55]. Thus, we choose the location of our “virtual” boundaries, l z , such that thecontributions of glial cell molecule uptake in y and spillover in z to the overall molecule concentrationdecay, λ , as defined in (23), are approximately equal and obtain l z = 2 µ m . According to [55,Fig. 5] this is indeed a biologically plausible value. Furthermore, to keep glial cell molecule uptakeconsistent with the default parameter values, the glial cell uptake rate constant is doubled to κ ,y = κ ,y = 5 . × − µ m µ s − .Since we are interested in the channel clearance due to glial cell uptake and there is a glial cellmembrane only in y , we consider the molecule uptake in y , obtained from (12) as i y ( t ) = ∞ (cid:88) n =1 ∞ (cid:88) m =1 ∞ (cid:88) p =1 A m,pn B m C p C ( x , y , z ) (cid:20) D sin( α m,pn l x ) + κ r α m,pn (1 − cos( α m,pn l x )) (cid:21) × [ κ ,y C m,y (0) + κ ,y C m,y ( l y )] (1 − cos( γ p l z )) γ p (1 − exp( − λ m,pn t )) λ m,pn . (55)as figure of interest.The results are depicted in Figs. 6 and 7.In Fig. 6, we consider different release positions by varying the release position in z , z in µ m . Notethat, in contrast to the results presented in Subsections VI-B and VI-C, the match between simulationdata and analytical results is not perfect. This is a consequence of the simulator design which doesnot allow for exact simulation of perfect adsorption, see [57] for a detailed analysis. However, theagreement between simulation data and analytical results is still very good given that the evaluatedquantities are cumulative, i.e., errors accumulate.In Fig. 7, the number of particles taken up at the glial cell after is shown depending on therelease position in z , z , for different values of y . First, we observe from Fig. 7 that the numberof particles taken up at the glial cell is independent of y . This is due to the fact that the reactionkinetics at the glial cell membrane are very slow compared to the diffusive motion of the particles; β is in fact very close to the reaction-induced lower bound from (48), such that the spatial concentrationprofile in y approaches the uniform distribution before the uptake exhibits a significant impact on themolecule concentration. Therefore, it is not relevant for the molecule uptake at the glial cell whetherparticles are released close to the glial cell membrane or not.In z direction, on the other hand, due to the instantaneously absorbing boundary, particle diffusionrather than reaction kinetics limits the number of escaping particles and thus also the number ofuptaken particles. Indeed, γ reaches the upper bound in (49) and hence the system is diffusion- N u m b e r o f U p t a k e n M o l ec u l e s i n y i y ( t ) [ - ] z = 1 . z = 1 . z = 1 . z = 1 . Figure 6. Number of uptaken particles at the glial cell in thepresence of spillover for different release positions z in µ m . Thesimulation data is depicted by diamond markers. z [ µ m] P a r t i c l e s U p t a k e n a t G li a l C e ll a f t e r m s i y ( ) [ - ] y = { . , . , . , . } Figure 7. Number of uptaken particles at the glial cell after in the presence of spillover for different release positions y , z in µ m . limited in z . In accordance with that, we observe from Figs. 6 and 7 that the number of particlestaken up at the glial cell heavily depends on the release position relative to the (virtual) systemboundaries in z . In particular, we observe from Fig. 7 that, as we increase or decrease z towards l z or , respectively, i y ( t ) decays the faster the closer the release site is to one of the boundaries,suggesting a nonlinear dependence of i y ( t ) on z .VII. C ONCLUSION
In this paper, we have proposed an analytical time domain expression for the CIR and the numberof uptaken particles in the tripartite synapse. The proposed model incorporates reversible bindingto individual postsynaptic receptors, molecule uptake at the presynaptic cell and at glial cells, andspillover to the extrasynaptic space. We have demonstrated the usefulness of our model by derivingfundamental bounds on the decay rate of the CIR in dependence on channel geometry and chemicalreaction kinetics. Furthermore, it was shown that the model allows novel quantitative insights regardingthe mitigation of ISI and the effectiveness of EH in synaptic communication assisted by a glial cell.Although we have presented some exemplary evaluations, we are convinced that, by informed use ofthis model, new input to many open questions regarding signal transmission in the tripartite synapsecan be generated and further contributions towards the design of synthetic neuronal communicationsystems and DMC systems in general are facilitated.Finally, we mention some limitations of the presented approach, which may motivate future research.First, the assumption that synaptic signal transmission can be modeled as a linear time-invariant systembreaks down as soon as any non-linear or time-dependent effects are considered. Thus, the model is not capable of describing short and long term synaptic plasticity or time-dependent release due toactivation of presynaptic auto-receptors. Second, the model that we used for postsynaptic receptors,is quite simplistic. A more elaborate model would consider receptors with multiple binding sites,multiple states, and state-dependent kinetics. Finally, it would be interesting to investigate the impactof receptor occupancy, i.e., the fact that a molecule cannot bind to a receptor binding side which isalready occupied by some other molecule. A PPENDIX
A. Derivation of the CIR in the Laplace Domain
Let us first consider the eigenfunctions { φ m , ψ p } of the Sturm-Liouville problems [58] ∂ ∂y φ = βφ, ≤ y ≤ l y , (56) ∂ ∂z ψ = γψ, ≤ z ≤ l z , (57)subject to k ,y ∂φ∂y = κ ,y φ at y = 0 k ,y ∂φ∂y = − κ ,y φ at y = l y , (58) k ,z ∂ψ∂z = κ ,z ψ at z = 0 k ,z ∂ψ∂z = − κ ,z ψ at z = l z , (59)where φ = φ ( y ) and ψ = ψ ( z ) are not identical . From [58, Thm. 3.3], we know that { φ m } and { ψ p } are countable sets and form an orthogonal basis of L ([0; l y ]) and L ([0; l z ]) , respectively, where L ([ a ; b ]) denotes the set of square-integrable functions on the interval [ a ; b ] , i.e., functions f , suchthat (cid:107) f (cid:107) a ; b ] := (cid:90) ba | f ( x ) | d x < ∞ . (60)Therefore, in particular, we can expand C in these bases. Assuming without loss of generality that { φ m } and { ψ p } are normalized, i.e., (cid:107) φ m (cid:107) y = 1 and (cid:107) ψ p (cid:107) z = 1 , this is done in the usual way byorthogonal projection [58], which is defined as P y g = ∞ (cid:88) m =1 b m φ m , P z h = ∞ (cid:88) p =1 c p ψ p , (61)where g ∈ L ([0; l y ]) , h ∈ L ([0; l z ]) , and the coefficients b m , c p are obtained as b m = (cid:90) l y gφ m d y, c p = (cid:90) l z hψ p d z. (62) Thus, if we apply (61) to (2) and exploit (56), (57), we obtain ∞ (cid:88) m =1 ∞ (cid:88) p =1 (cid:20) ∂ C m,p ∂x − β m C m,p − γ p C m,p − D ∂C m,p ∂t (cid:21) φ m ψ p = 0 , (63)where the coefficients C m,p = C m,p ( x, t ) depend on x and t . Now, because the { φ m } are linearlyindependent and also the { ψ p } are linearly independent, (63) implies that ∂ C m,p ∂x − β m C m,p − γ p C m,p − D ∂C m,p ∂t = 0 ∀ m, p. (64)Furthermore, C m,p is to fulfill boundary conditions (3), (4) and initial condition C m,p, = C m,p ( x,
0) = (cid:90) l z (cid:90) l y C φ m ψ p d y d z = δ ( x − x ) φ m ( y ) ψ p ( z ) . (65)Now, we apply the Laplace transform with respect to t , denoted by L , to (64) and obtain theinhomogeneous ordinary differential equation ∂ ¯ C m,p ∂x = (cid:16) sD + β m + γ p (cid:17) ¯ C m,p − D C m,p, , (66)where s denotes the Laplace variable, ¯ C m,p ( x, s ) = L{ C m,p ( x, t ) } = (cid:82) ∞ C m,p ( x, τ ) exp( − sτ )d τ , and(66) is subject to D ∂ ¯ C m,p ∂x = κ r ¯ C m,p , at x = 0 , and (67) D ∂ ¯ C m,p ∂x = − κ a ¯ C m,p − κ d D s ∂ ¯ C m,p ∂x , at x = l x . (68)Here, we have exploited the properties of the Laplace transform that L{ d f d t } = s L{ f } − f (0) and L{ (cid:82) t f d τ } = s L{ f } . To solve (66)–(68), we follow a similar approach as [40, Ch. 14] and decompose C m,p as C m,p = U m,p + W m,p , (69)such that both, U m,p and W m,p , fulfill (64), U m,p = C m,p, at t = 0 , W m,p = 0 at t = 0 , and the sumof U m,p and W m,p fulfills (3) and (4). Let us denote the Laplace transforms of U m,p and W m,p as ¯ U m,p and ¯ W m,p , respectively, and define q (cid:48) = (cid:112) sD + β m + γ p . Then, from [40, Ch. 14.3], we find that ¯ U m,p = φ m ( y ) ψ p ( z )2 Dq (cid:48) e − q (cid:48) | x − x | , (70) ¯ C m,p = ¯ U m,p + ¯ W m,p = φ m ( y ) ψ p ( z )2 Dq (cid:48) e − q (cid:48) | x − x | + A sinh( q (cid:48) x ) + B cosh( q (cid:48) x ) , (71) where A and B are constants with respect to x and must be chosen such that (67) and (68) arefulfilled. Substituting (71) in (67) and (68) yields A = 1 q (cid:48) (cid:18) K ,x ¯ U m,p (cid:12)(cid:12) x =0 − ∂ ¯ U m,p ∂x (cid:12)(cid:12)(cid:12)(cid:12) x =0 + BK ,x (cid:19) , (72) B = (cid:20) − q (cid:48) ∂ ¯ U m,p ∂x (cid:12)(cid:12)(cid:12)(cid:12) x = l x + q (cid:48) (cid:18) ∂ ¯ U m,p ∂x (cid:12)(cid:12)(cid:12)(cid:12) x =0 − K ,x ¯ U m,p (cid:12)(cid:12) x =0 (cid:19) cosh( q (cid:48) l x ) − q (cid:48) K ,x ¯ U m,p (cid:12)(cid:12) x = l x − (cid:18) K ,x K ,x ¯ U m,p (cid:12)(cid:12) x =0 − K ,x ∂ ¯ U m,p ∂x (cid:12)(cid:12)(cid:12)(cid:12) x =0 (cid:19) sinh( q (cid:48) l x ) (cid:21) / (cid:2) q (cid:48) ( K ,x + K ,x ) cosh( q (cid:48) l x ) + (cid:0) ( q (cid:48) ) + K ,x K ,x (cid:1) sinh( q (cid:48) l x ) (cid:3) , (73)and we have set q (cid:48) = (cid:112) sD + β m + γ p , K ,x = κ r /D and K ,x = κ a / (cid:2) D (1 + κ d s ) (cid:3) . B. Time Domain Solution
Now, to transform (71) back to time domain, we need to compute the inverse Laplace transform, C m,p ( x, t ) = 12 πj (cid:90) ν + j ∞ ν − j ∞ e st ¯ C m,p ( x, s ) d s, (74)where j denotes the imaginary unit and ν ∈ R + needs to be chosen large enough such that allsingularities of ¯ C m,p are left of the line along which the integral is computed.Similarly to [40], instead of integrating along the infinite line, we complete the integration pathwith a semicircle to a closed contour which contains the origin and all singularities of ¯ C m,p . Then,we use the residue theorem to replace this contour integral with a sum over the residues { σ } of ¯ C m,p , C m,p ( x, t ) = (cid:88) { σ } Res( e st ¯ C m,p , σ ) . (75)Now, (71) has poles at q (cid:48) = 0 and infinitely many non-zero simple poles defined as the solutionsof ( K ,x + K ,x ) cosh( q (cid:48) l x ) + (cid:18) ( q (cid:48) ) + K ,x K ,x q (cid:48) (cid:19) sinh( q (cid:48) l x ) = 0 . (76)Note that, at q (cid:48) = 0 , (76) does not have a pole but a removable singularity as lim q (cid:48) → sinh( q (cid:48) l x ) /q (cid:48) = 1 .If K ,x = 0 and β m = γ p = 0 , (76) has a root at q (cid:48) = 0 and thus contributes a pole at q (cid:48) = 0 to (71).To compute (75), we start with the non-zero roots of (76).The residues of e st ¯ C m,p at simple poles σ can be computed with l’Hôpital’s rule asRes (cid:0) e st ¯ C m,p , σ (cid:1) = e st ¯ C Nm,p∂ ¯ C Dm,p ∂s (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = σ , (77)where ¯ C Nm,p and ¯ C Dm,p denote the numerator and denominator of ¯ C m,p , respectively, as defined by (71). Now, we set q (cid:48) = jα and from (76) it follows that the non-zero poles of (71) correspond to thenon-zero roots α m,pn , n ∈ N of (24) with either positive real part or positive imaginary part. Eq. (24)admits infinitely many pairs of real solutions and, under some conditions (if κ d > and β m + γ p is large enough), also one pair of complex conjugate pure imaginary roots. Due to our substitution q (cid:48) = (cid:112) sD + β m + γ p , we need to consider only one root from each pair of roots and w.l.o.g. wechoose the one with positive real or positive imaginary part, respectively.Furthermore, we index the real roots ascendingly, starting with index if (24) does not admitimaginary roots and starting with index otherwise. In the latter case, we index the imaginary rootof (24) with index .Applying (77) to (71) we obtainRes (cid:0) e st ¯ C m,p , σ m,pn (cid:1) = φ m ( y ) ψ p ( z )( κ r sin( α m,pn x ) + Dα m,pn cos( α m,pn x )) (78) × α m,pn ) D ( λ m,pn − κ d ) + ( λ m,pn ) κ a ] / { ( λ m,pn − κ d )[ D ( α m,pn ) ( κ r ( λ m,pn − κ d ) + κ a λ m,pn ) + κ r κ a Dλ m,pn ]++ 2 κ r κ a κ d D ( α m,pn ) + κ r κ a D ( λ m,pn ) + 2 κ a κ d D ( α m,pn ) ++ l x [ κ r D ( α m,pn ) ( λ m,pn − κ d ) + κ a D ( α m,pn ) ( λ m,pn ) ++ D ( α m,pn ) ( λ m,pn − κ d ) + κ r κ a ( λ m,pn ) ] } (79) × ( κ r sin( α m,pn x ) + Dα m,pn cos( α m,pn x )) exp( − λ m,pn t ) , (80)where σ m,pn = − (( α m,pn ) + β m + γ p ) D and λ m,pn is defined in (23).Now, for the poles of (71) at q (cid:48) = 0 , we employ a different strategy. Namely, here we use the factthat the residue of a function at a pole σ is given as the coefficient of the term of order − of itsLaurent series expansion at σ . Hence, if the Laurent series expansion of e st ¯ C m,p at σ is given as e st ¯ C m,p = ∞ (cid:88) k = −∞ c k ( s − σ ) k , (81)we have Res ( e st ¯ C m,p , σ ) = c − .For q (cid:48) = 0 , we have s = − ( β m + γ p ) D and thus we need the Laurent series expansion of e st ¯ C m,p at σ = − ( β m + γ p ) D . Before we perform the series expansion, we apply the change of variable q (cid:48) = (cid:112) sD + β m + γ p and note that ( s − σ ) − = ( s + ( β m + γ p ) D ) − = D − ( q (cid:48) ) − . (82)Thus, we can equivalently find the series expansion of e st ¯ C m,p at q (cid:48) = 0 , e st ¯ C m,p = ∞ (cid:88) k = −∞ d k ( q (cid:48) ) k , (83) and have Res ( e st ¯ C m,p , σ ) = c − = d − D .First, if K ,x (cid:54) = 0 , (71) has a simple pole at q (cid:48) = 0 and by expanding the expression in q (cid:48) = 0 , itcan easily be shown that in this case Res ( e st ¯ C m,p , σ ) = 0 . The same observation holds, if K ,x = 0 but β m + γ p (cid:54) = 0 .Now, if K ,x = 0 and β m + γ p = 0 , there are non-zero contributions from only two terms in (71),namely, ∂ ¯ U m,p ∂x (cid:12)(cid:12)(cid:12) x =0 cosh( q (cid:48) l x ) cosh( q (cid:48) x ) q (cid:48) ( K ,x q (cid:48) cosh( q (cid:48) l x ) + sinh( q (cid:48) l x )) (84)and − ∂ ¯ U m,p ∂x (cid:12)(cid:12)(cid:12) x = l x cosh( q (cid:48) x ) q (cid:48) ( K ,x q (cid:48) cosh( q (cid:48) l x ) + sinh( q (cid:48) l x )) . (85)For (84), with s = ( q (cid:48) ) D , we have ∂ ¯ U m,p ∂x (cid:12)(cid:12)(cid:12) x =0 cosh( q (cid:48) l x ) cosh( q (cid:48) x ) e ( q (cid:48) ) Dt q (cid:48) ( K ,x q (cid:48) cosh( q (cid:48) l x ) + sinh( q (cid:48) l x )) = φ m ( y ) ψ p ( z )2 D ( q (cid:48) ) cosh( q (cid:48) l x ) cosh( q (cid:48) x ) e ( q (cid:48) ) Dt − q (cid:48) x ( K ,x ( q (cid:48) ) cosh( q (cid:48) l x ) + sinh( q (cid:48) l x ) q (cid:48) ) . (86)Now, we denote the Taylor series expansion of the right-hand side of (86) at q (cid:48) = 0 as cosh( q (cid:48) l x ) cosh( q (cid:48) x ) e q (cid:48) ( q (cid:48) Dt − x ) K ,x ( q (cid:48) ) cosh( q (cid:48) l x ) + sinh( q (cid:48) l x ) q (cid:48) = d (1)0 + d (1)1 q (cid:48) + d (1)2 ( q (cid:48) ) + . . . (87)and find O ( q (cid:48) ) = cosh( q (cid:48) l x ) cosh( q (cid:48) x ) e q (cid:48) ( q (cid:48) Dt − x ) = ( d (1)0 + O ( q (cid:48) )) (cid:18) K ,x ( q (cid:48) ) cosh( q (cid:48) l x ) + sinh( q (cid:48) l x ) q (cid:48) (cid:19) = ( d (1)0 + O ( q (cid:48) )) (cid:40) κ a κ d ( q (cid:48) ) Dκ d + 1 (cid:2) O (( q (cid:48) ) ) (cid:3) + (cid:2) l x + O (( q (cid:48) ) ) (cid:3)(cid:41) = ( d (1)0 + O ( q (cid:48) )) (cid:26) κ a κ d (cid:2) O (( q (cid:48) ) ) (cid:3) (cid:2) O (( q (cid:48) ) ) (cid:3) + (cid:2) l x + O (( q (cid:48) ) ) (cid:3)(cid:27) , (88)where O ( · ) denotes big O notation and we have used the well-known series expansions of e x , sinh( x ) , cosh( x ) , and / ( x + 1) . Collecting the constant terms in (88), we find d (1)0 (cid:18) κ a κ d + l x (cid:19) , (89)and therefore d (1)0 = κ d κ a + l x κ d . (90) Multiplying both sides of (88) by φ m ( y ) ψ p ( z ) / (2 D ( q (cid:48) ) ) , we obtain the Laurent series expansionof the left-hand side of (86) as ∂ ¯ U m,p ∂x (cid:12)(cid:12)(cid:12) x =0 cosh( q (cid:48) l x ) cosh( q (cid:48) x ) e ( q (cid:48) ) Dt q (cid:48) ( K ,x q (cid:48) cosh( q (cid:48) l x ) + sinh( q (cid:48) l x )) = κ d κ a + l x κ d φ m ( y ) ψ p ( z )2 D ( q (cid:48) ) + O (( q (cid:48) ) − ) (91)and in a similar fashion we obtain − ∂ ¯ U m,p ∂x (cid:12)(cid:12)(cid:12) x = l x cosh( q (cid:48) x ) e ( q (cid:48) ) Dt q (cid:48) ( K ,x q (cid:48) cosh( q (cid:48) l x ) + sinh( q (cid:48) l x )) = κ d κ a + l x κ d φ m ( y ) ψ p ( z )2 D ( q (cid:48) ) + O (( q (cid:48) ) − ) . (92)Finally, summing (91) and (92), we conclude that d − = κ d ( φ m ( y ) ψ p ( z )) / [ D ( κ a + l x κ d )] and thus c − = κ d ( φ m ( y ) ψ p ( z )) / ( κ a + l x κ d ) .Putting together the pieces, we obtain C m,p ( x, t ) = κ d κ a + l x κ d φ m ( y ) ψ p ( z ) { } ( κ r ) { } ( β m ) { } ( γ p ) (93) + ∞ (cid:88) n =1 A m,pn C ( x , y , z ) C m,pn,x ( x ) exp( − λ m,pn t ) , (94)where A m,pn , C ( x , y , z ) , and C m,pn,x ( x ) are defined in (20), (19), and (16), respectively.Now, to obtain the final solution C ( x, y, z, t ) as a linear combination of orthonormal basis functions,we need to determine { φ m } and { ψ p } . But these are well-known in the literature, see for example[40], [58], and we repeat them here for completeness only: φ m ( y ) = (cid:40) (cid:2) − { } ( β m ) (cid:3) ( β m k ,y + κ ,y ) (cid:2) l y ( κ ,y + k ,y β m ) + κ ,y k ,y (cid:3) ( κ ,y + k ,y β m ) + κ ,y k ,y ( κ ,y + k ,y β m ) (cid:41) × [ κ ,y sin( β m y ) + k ,y β m cos( β m y )] , (95) ψ p ( z ) = (cid:40) (cid:2) − { } ( γ p ) (cid:3) ( γ p k ,z + κ ,z ) (cid:2) l z ( κ ,z + k ,z γ p ) + κ ,z k ,z (cid:3) ( κ ,z + k ,z γ p ) + κ ,z k ,z ( κ ,z + k ,z γ p ) (cid:41) × [ κ ,z sin( γ p z ) + k ,z γ p cos( γ p z )] , (96)where β m and γ p are the positive, or, if κ ,y = κ ,y = 0 ( κ ,z = κ ,z = 0) , the non-negative, rootsof (25), (26). All special cases, e.g. perfectly absorbing boundaries k ,y = k ,y = k ,z = k ,z = 0 ,follow from this definition. In some special cases, (25), (26) simplify considerably and an explicitexpressions can be given for { β m } and { γ p } . For a comprehensive listing of different combinations ofboundary conditions and the associated eigenvalues, see for example [38, Tab. 2.1]. If the boundary in y ( z ) is reflective, i.e., κ ,y = κ ,y = 0 ( κ ,z = κ ,z = 0) , the constant eigenfunction / (cid:112) l y ( / √ l z ) isadmitted by (56), (58) ((57), (59)), and thus also the eigenvalues β = 0 ( γ = 0 ) need to be considered.This completes the proof. R EFERENCES [1] S. Lotter, A. Ahmadzadeh, and R. Schober, “Channel modeling for synaptic molecular communication with re-uptake and reversiblereceptor binding,” accepted for presentation at IEEE Int. Conf. Commun. (ICC), 2020.[2] T. Nakano, A. W. Eckford, and T. Haraguchi,
Molecular Communication . Cambridge University Press, 2013.[3] I. F. Akyildiz, M. Pierobon, S. Balasubramaniam, and Y. Koucheryavy, “The internet of bio-nano things,”
IEEE Commun. Mag. ,vol. 53, no. 3, pp. 32–40, Mar. 2015.[4] R. A. Freitas,
Nanomedicine, Volume I: Basic Capabilities . Landes Bioscience Georgetown, TX, 1999, vol. 1.[5] M. Veleti´c and I. Balasingham, “Synaptic communication engineering for future cognitive brain–machine interfaces,”
Proc. IEEE ,vol. 107, no. 7, pp. 1425–1441, Jul. 2019.[6] D. W. K. Ng, E. S. Lo, and R. Schober, “Wireless information and power transfer: Energy efficiency optimization in OFDMAsystems,”
IEEE Trans. Wireless Commun. , vol. 12, no. 12, pp. 6352–6370, Dec. 2013.[7] Y. Deng, W. Guo, A. Noel, A. Nallanathan, and M. Elkashlan, “Enabling energy efficient molecular communication via moleculeenergy transfer,”
IEEE Commun. Lett. , vol. 21, no. 2, pp. 254–257, Feb. 2017.[8] W. Guo, Y. Deng, H. B. Yilmaz, N. Farsad, M. Elkashlan, A. Eckford, A. Nallanathan, and C. Chae, “SMIET: Simultaneousmolecular information and energy transfer,”
IEEE Wireless Commun. , vol. 25, no. 1, pp. 106–113, Feb. 2018.[9] V. Musa, G. Piro, L. A. Grieco, and G. Boggia, “A lean control theoretic approach to energy-harvesting in diffusion-basedmolecular communications,”
IEEE Commun. Lett. , pp. 1–1, Feb. 2020.[10] H. Arjmandi, M. Movahednasab, A. Gohari, M. Mirmohseni, M. Nasiri-Kenari, and F. Fekri, “ISI-avoiding modulation fordiffusion-based molecular communication,”
IEEE Trans. Mol. Biol. Multi-Scale Commun. , vol. 3, no. 1, pp. 48–59, Mar. 2017.[11] M. S. Leeson and M. D. Higgins, “Forward error correction for molecular communications,”
Nano Commun. Networks , vol. 3,no. 3, pp. 161–167, Sep. 2012.[12] P. Shih, C. Lee, P. Yeh, and K. Chen, “Channel codes for reliability enhancement in molecular communication,”
IEEE J. Sel.Areas Commun. , vol. 31, no. 12, pp. 857–867, Dec. 2013.[13] B. Tepekule, A. E. Pusane, H. B. Yilmaz, C. Chae, and T. Tugcu, “ISI mitigation techniques in molecular communication,”
IEEETrans. Mol. Biol. Multi-Scale Commun. , vol. 1, no. 2, pp. 202–216, Jun. 2015.[14] A. Noel, K. C. Cheung, and R. Schober, “Improving receiver performance of diffusive molecular communication with enzymes,”
IEEE Trans. Nanobiosci. , vol. 13, no. 1, pp. 31–43, Mar. 2014.[15] A. C. Heren, H. B. Yilmaz, C. Chae, and T. Tugcu, “Effect of degradation in molecular communication: Impairment orenhancement?”
IEEE Trans. Mol. Biol. Multi-Scale Commun. , vol. 1, no. 2, pp. 217–229, Jun. 2015.[16] L. K. Bak, A. Schousboe, and H. S. Waagepetersen, “The glutamate/GABA-glutamine cycle: Aspects of transport, neurotransmitterhomeostasis and ammonia transfer,”
J. Neurochem. , vol. 98, no. 3, pp. 641–653, Jun. 2006.[17] T. A. Nielsen, D. A. DiGregorio, and R. Silver, “Modulation of glutamate mobility reveals the mechanism underlying slow-risingAMPAR EPSCs and the diffusion coefficient in the synaptic cleft,”
Neuron , vol. 42, no. 5, pp. 757 – 771, Jun. 2004.[18] J. Hasegawa, T. Obara, K. Tanaka, and M. Tachibana, “High-density presynaptic transporters are required for glutamate removalfrom the first visual synapse,”
Neuron , vol. 50, no. 1, pp. 63–74, Apr. 2006.[19] N. C. Danbolt, “Glutamate uptake,”
Prog. Neurobiol. , vol. 65, no. 1, pp. 1–105, Sep. 2001.[20] N. J. Maragakis and J. D. Rothstein, “Glutamate transporters in neurologic disease,”
Arch. Neurol. , vol. 58, no. 3, pp. 365–370,Mar. 2001.[21] L. Hertz, R. Dringen, A. Schousboe, and S. R. Robinson, “Astrocytes: Glutamate producers for neurons,”
J. Neurosci. Res. ,vol. 57, no. 4, pp. 417–428, Oct. 1999.[22] E. Balevi and O. B. Akan, “A physical channel model for nanoscale neuro-spike communications,”
IEEE Trans. Commun. , vol. 61,no. 3, pp. 1178–1187, Mar. 2013. [23] M. Veleti´c, F. Mesiti, P. A. Floor, and I. Balasingham, “Communication theory aspects of synaptic transmission,” in Proc. IEEEIntern. Conf. Commun. , 2015, pp. 1116–1121.[24] T. Khan, B. A. Bilgin, and O. B. Akan, “Diffusion-based model for synaptic molecular communication channel,”
IEEE Trans.Nanobiosci. , vol. 16, no. 4, pp. 299–308, Jun. 2017.[25] S. Nadkarni and P. Jung, “Modeling synaptic transmission of the tripartite synapse,”
Phys. Biol. , vol. 4, no. 1, pp. 1–9, Jan. 2007.[26] A. Mahdavi, F. Bahrami, and M. Janahmadi, “Analysis of impaired LTP in schizophrenia using an extended mathematical modelof a tripartite synapse,” in , Nov. 2015, pp. 76–80.[27] F. Mesiti, P. A. Floor, and I. Balasingham, “Astrocyte to neuron communication channels with applications,”
IEEE Trans. Mol.Biol. Multi-Scale Commun. , vol. 1, no. 2, pp. 164–175, Jun. 2015.[28] M. G. Stewart, V. I. Popov, I. V. Kraev, N. Medvedev, and H. A. Davies, “Chapter one - structure and complexity of the synapseand dendritic spine,” in
The Synapse , V. Pickel and M. Segal, Eds. Boston: Academic Press, 2014, pp. 1 – 20.[29] A. Noel, “Modeling and analysis of diffusive molecular communication systems,” Ph.D. dissertation, University of BritishColumbia, 2015.[30] E. L. Cussler,
Diffusion: Mass Transfer in Fluid Systems . Cambridge University Press, 2009.[31] A. Ahmadzadeh, H. Arjmandi, A. Burkovski, and R. Schober, “Comprehensive reactive receiver modeling for diffusive molecularcommunication systems: Reversible binding, molecule degradation, and finite number of receptors,”
IEEE Trans. Nanobiosci. ,vol. 15, no. 7, pp. 713–727, Oct. 2016.[32] M. Sheng and C. C. Hoogenraad, “The postsynaptic architecture of excitatory synapses: A more quantitative view,”
Annu. Rev.Biochem. , vol. 76, no. 1, pp. 823–847, 2007, pMID: 17243894.[33] R. Zwanzig and A. Szabo, “Time dependent rate of diffusion-influenced ligand binding to receptors on cell surfaces,”
Biophys.J. , vol. 60, no. 3, pp. 671–678, Sep. 1991.[34] A. M. Berezhkovskii, Y. A. Makhnovskii, M. I. Monine, V. Y. Zitserman, and S. Y. Shvartsman, “Boundary homogenization fortrapping by patchy surfaces,”
J. Chem. Phys. , vol. 121, no. 22, pp. 11 390–11 394, Nov. 2004.[35] A. Szabo, “Theoretical approaches to reversible diffusion-influenced reactions: Monomer–excimer kinetics,”
J. Chem. Phys. ,vol. 95, no. 4, pp. 2481–2490, May 1991.[36] V. I. Zarnitsyna, J. Huang, F. Zhang, Y.-H. Chien, D. Leckband, and C. Zhu, “Memory in receptor–ligand-mediated cell adhesion,”
Proc. Nat. Acad. Sci. , vol. 104, no. 46, pp. 18 037–18 042, Nov. 2007.[37] F. C. Collins and G. E. Kimball, “Diffusion-controlled reaction rates,”
J. Colloid Sci. , vol. 4, no. 4, pp. 425–437, Aug. 1949.[38] M. Ozisik,
Boundary Value Problems of Heat Conduction , ser. Dover Books on Engineering. Dover Publications, 2013.[39] N. Agmon, “Diffusion with back reaction,”
J. Chem. Phys. , vol. 81, no. 6, pp. 2811–2817, Apr. 1984.[40] H. Carslaw and J. Jaeger,
Conduction of Heat in Solids , ser. Oxford Science Publications. Clarendon Press, 1986.[41] R. Heidelberger, W. B. Thoreson, and P. Witkovsky, “Synaptic transmission at retinal ribbon synapses,”
Prog. Retin. Eye Res. ,vol. 24, no. 6, pp. 682–720, Nov. 2005.[42] S. S. Andrews, “Accurate particle-based simulation of adsorption, desorption and partial transmission,”
Phys. Biol. , vol. 6, no. 4,p. 046015, Nov. 2009.[43] H. Berg and E. Purcell, “Physics of chemoreception,”
Biophys. J. , vol. 20, no. 2, pp. 193 – 219, Nov. 1977.[44] D. Shoup and A. Szabo, “Role of diffusion in ligand binding to macromolecules and cell-bound receptors,”
Biophys. J. , vol. 40,no. 1, pp. 33–39, Oct. 1982.[45] I. V. Gopich and A. Szabo, “Kinetics of reversible diffusion influenced reactions: The self-consistent relaxation time approxima-tion,”
J. Chem. Phys. , vol. 117, no. 2, pp. 507–517, Apr. 2002.[46] A. M. Berezhkovskii and A. Szabo, “Effect of ligand diffusion on occupancy fluctuations of cell-surface receptors,”
J. Chem.Phys. , vol. 139, no. 12, p. 121910, Jul. 2013. [47] A. M. Berezhkovskii, L. Dagdug, V. A. Lizunov, J. Zimmerberg, and S. M. Bezrukov, “Communication: Clusters of absorbingdisks on a reflecting wall: Competition for diffusing particles,” J. Chem. Phys. , vol. 136, no. 21, p. 211102, May 2012.[48] A. Mugler, A. G. Bailey, K. Takahashi, and P. R. ten Wolde, “Membrane clustering and the role of rebinding in biochemicalsignaling,”
Biophys. J. , vol. 102, no. 5, pp. 1069–1078, Mar. 2012.[49] U. Ashery, N. Bielopolski, A. Lavi, B. Barak, L. Michaeli, Y. Ben-Simon, A. Sheinin, D. Bar-On, Z. Shapira, and I. Gottfried,“Chapter two - the molecular mechanisms underlying synaptic transmission: A view of the presynaptic terminal,” in
The Synapse ,V. Pickel and M. Segal, Eds. Boston: Academic Press, 2014, pp. 21 – 109.[50] W. Holmes, “Modeling the effect of glutamate diffusion and uptake on NMDA and non-NMDA receptor saturation,”
Biophys. J. ,vol. 69, no. 5, pp. 1734 – 1747, Nov. 1995.[51] P. Jonas, G. Major, and B. Sakmann, “Quantal components of unitary EPSCs at the mossy fibre synapse on CA3 pyramidal cellsof rat hippocampus.”
J. Physiol. (Lond.) , vol. 472, no. 1, pp. 615–663, Dec. 1993.[52] J. I. Wadiche and M. P. Kavanaugh, “Macroscopic and microscopic properties of a cloned glutamate transporter/chloride channel,”
J. Neurosci. , vol. 18, no. 19, pp. 7650–7661, Oct. 1998.[53] K. P. Lehre and N. C. Danbolt, “The number of glutamate transporter subtype molecules at glutamatergic synapses: Chemicaland stereological quantification in young adult rat brain,”
J. Neurosci. , vol. 18, no. 21, pp. 8751–8757, Nov. 1998.[54] L. P. Savtchenko, S. Sylantyev, and D. A. Rusakov, “Central synapses release a resource-efficient amount of glutamate,”
Nat.Neurosci. , vol. 16, no. 1, p. 10, Dec. 2012.[55] R. Ventura and K. M. Harris, “Three-dimensional relationships between hippocampal synapses and astrocytes,”
J. Neurosci. ,vol. 19, no. 16, pp. 6897–6906, Aug. 1999.[56] A. M. Sweeney, K. E. Fleming, J. P. McCauley, M. F. Rodriguez, E. T. Martin, A. A. Sousa, R. D. Leapman, and A. Scimemi,“PAR1 activation induces rapid changes in glutamate uptake and astrocyte morphology,”
Sci. Rep. , vol. 7, p. 43606, Jan. 2017.[57] P. Clifford and N. Green, “On the simulation of the Smoluchowski boundary condition and the interpolation of Brownian paths,”
Molecular Physics , vol. 57, no. 1, pp. 123–128, Jan. 1986.[58] G. Cain and G. H. Meyer,