Channel Modeling for Synaptic Molecular Communication With Re-uptake and Reversible Receptor Binding
CChannel Modeling for Synaptic MolecularCommunication With Re-uptake and ReversibleReceptor Binding
Sebastian Lotter, Arman Ahmadzadeh, and Robert Schober
Friedrich-Alexander University Erlangen-Nuremberg, Germany
Abstract —In Diffusive Molecular Communication (DMC), in-formation is transmitted by diffusing molecules. Synaptic signal-ing is a natural implementation of this paradigm. It is responsiblefor relaying information from one neuron to another, but alsoprovides support for complex functionalities, such as learningand memory. Many of its features are not yet understood,some are, however, known to be critical for robust, reliableneural communication. In particular, some synapses feature are-uptake mechanism at the presynaptic neuron, which providesa means for removing neurotransmitters from the synapticcleft and for recycling them for future reuse. In this paper,we develop a comprehensive channel model for synaptic DMCencompassing a spatial model of the synaptic cleft, moleculere-uptake at the presynaptic neuron, and reversible binding toindividual receptors at the postsynaptic neuron. Based on thismodel, we derive an analytical time domain expression for thechannel impulse response (CIR) of the synaptic DMC system.Our model explicitly incorporates macroscopic physical channelparameters and can be used to evaluate the impact of re-uptake,receptor density, and channel width on the CIR of the synapticDMC system. Furthermore, we provide results from particle-based computer simulation, which validate the analytical model.The proposed comprehensive channel model for synaptic DMCsystems can be exploited for the investigation of challengingproblems, like the quantification of the inter-symbol interferencebetween successive synaptic signals and the design of syntheticneural communication systems.
I. I
NTRODUCTION
In nature, information exchange between many differentbiological entities, such as cells, organs, or even differentindividuals, is based on the release, propagation, and sensingof molecules. This process is called Molecular Communica-tion (MC). Although traditionally studied by biologists andmedical scientists, it has recently also attracted interest in thecommunications research community [1]. MC is envisioned toopen several exciting new application areas for which commu-nication at nano-scale is key [2]. One of these applications isthe deployment and operation of synthetic cells in the humanbody for the purpose of tumor detection or treatment [3],another one is the development of brain-machine interfaces(BMIs) for detection or replacement of dysfunctional neuralunits [4]. In both cases, communication at cell level, eitherbetween synthetic cells or between synthetic and natural cells,is required, but can not be implemented using traditionalwireless communication systems.MC systems in which molecules propagate via Brownianmotion, referred to as Diffusive Molecular Communication(DMC) systems, provide a viable alternative, as neither specialinfrastructure, nor external energy supply is needed. The design of synthetic DMC systems, however, poses severalchallenges. Firstly, in environments for which synthetic DMCis intended, such as the human body, the amount of moleculesavailable for transmission is typically limited. Thus, the em-ployed communication scheme needs to be extremely energyefficient. The concept of Energy Harvesting (EH) enablesultra-low-power applications in traditional communications [5]and several methods to leverage it also for DMC have beenproposed recently [6], [7]. The models considered in [6], [7],however, are based on free-space propagation and thus do nottake into account a particular communication environment.Secondly, because Brownian motion is an undirected prop-agation mechanism, DMC channels are typically dispersiveand adequate measures to mitigate inter-symbol interference(ISI) need to be taken. Several approaches for ISI mitigationin DMC have been proposed in the literature, including ISI-aware modulation schemes [8], forward error-correction codes[9], [10], and equalization techniques [11].A difficulty in resolving these issues is that existing conceptsfrom traditional communications cannot be easily transferredto MC due to limited processing capabilities at both transmitterand receiver. On the other hand, natural DMC systems haveevolved over millions of years to deal with such challenges.Inspired by such natural systems, enzymatic degradation ofinformation molecules in the channel is considered for themitigation of ISI in [12], [13]. This approach is interesting,because it directly impacts the channel characteristics and doesnot increase the complexity of the transmitter or receiver. Itdoes, however, incur higher energy cost for the production ofenzymes and additional signaling molecules, and is thus notnecessarily energy efficient.While the approach in [12] was inspired by the neu-romuscular junction, we consider a different natural DMCsystem in this paper, namely molecular synaptic transmissionbetween two neurons. Abstracted in communication terms,here, the presynaptic neuron (transmitter) encodes a sequenceof electrical impulses (data stream) into a spatio-temporal neu-rotransmitter release pattern (molecular signal) which prop-agates through the synaptic cleft (channel) and is finallyreceived by the postsynaptic neuron (receiver) where theneurotransmitters activate membrane receptors, see Fig. 1. Inaddition, transporter proteins at the presynaptic neuron providea re-uptake mechanism for many common neurotransmitters,such as dopamine, serotonin, norepinephrine, γ -aminobutyricacid (GABA), and glycine [14]. In this way, the channel a r X i v : . [ q - b i o . S C ] D ec s cleared and signaling molecules are recycled for futurereuse. It is known that this re-uptake mechanism is criticalfor synaptic communication; several severe mental diseases,including attention deficit hyperactivity disorder and epilepsyare associated with dysfunctional re-uptake [14].Despite its importance for neural information transmissionand its bio-physical characteristics, most existing models ofsynaptic signaling are based on free-space propagation [15],[16] and can thus not capture the impact of specific synapticchannel parameters. Recent progress in this direction hasbeen reported in [17]. In this article, a new analytical modelfor molecule diffusion in the synaptic cleft, incorporating asimplified geometric representation of the synaptic cleft asinfinite region bounded by to parallel planes, is proposed. Themodel in [17] does also take into account molecule re-uptakeand binding to postsynaptic receptors, but the iterative schemeused to determine postsynaptic binding does not lend itselfto analytical investigations. A closed-form expression for thechannel impulse response (CIR) is not provided. Furthermore,to evaluate the synaptic CIR with respect to ISI, moleculerebinding at the postsynaptic neuron, i.e., reversible binding,needs to be considered, while [17] assumes irreversible bind-ing.To the best of the authors’ knowledge, there is no chan-nel model available for the synaptic cleft, which allows toassess the quantitative impact of molecule re-uptake and otherchannel parameters on the CIR. Such a model would allowthe investigation of the relevance of molecule re-uptake forEH and ISI mitigation, and on the other hand, be of interestin its own right, e.g. for the development of brain-machineinterfaces.The main contribution of this paper is an analytical timedomain expression for the CIR of the synaptic cleft whichincorporates biologically plausible models of presynaptic re-uptake, postsynaptic reversible binding kinetics, and the partic-ular cleft geometry. Furthermore, this expression is validatedby particle-based computer simulation and experimentallyshown to provide an unbiased estimator for the CIR.The remainder of this paper is organized as follows: In Sec-tion II, we state the system model and the main assumptionsused in Section III to derive the CIR of the synaptic cleft. Theparticle-based simulator design is outlined in Section IV, andnumerical results are presented in Section V. Finally, the mainfindings are summarized in Section VI.II. S YSTEM M ODEL
A. Geometry and Assumptions
The complex geometry of the synaptic cleft is abstractedin our model as a 3-dimensional rectangular cuboid, withfaces in x -direction representing the membranes of the pre-and postsynaptic neurons (see Fig. 1). This model is similarto the geometric model used in [17], with the differencethat the model considered here is bounded in all dimensions,while in [17], y and z extend to infinity. This abstraction isanalytically more tractable than the actual non-regular shapedsynaptic domain, while still retaining its characteristic features.In the context of this model, we will refer to the boundary representing the pre- and postsynaptic membranes as the left and right boundaries, respectively. Formally, we denote thedomain of the synaptic cleft in Cartesian coordinates as Ω = { ( 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 atany location within the box defined by (1) as C Ω ( x, y, z, t ) .To derive the (dimensionless) impulse response of thesynaptic channel, h ( t ) , we consider instantaneous release of N particles at time t = t at location ( x, y, z ) = ( x , y , z ) ∈ Ω .Without loss of generality (w.l.o.g.), we set t = 0 . h ( t ) is thengiven as the number of particles adsorbed to the receptors ofthe postsynaptic membrane as a function of t, t ≥ . For ouranalysis, we make four assumptions:A1) The faces in y and z direction are fully reflective.A2) The receptors at the postsynaptic membrane are uni-formly distributed.A3) The receptors cannot be occupied, i.e., multiple particlesmay bind to a single receptor.A4) Reversible adsorption to individual receptors with intrin-sic association coefficient κ a in µ m µ s − and intrinsicdissociation rate κ d in µ s − can be treated equivalentlyas reversible adsorption to a homogeneous surface witheffective association coefficient κ a in µ m µ s − and dis-sociation rate κ d .A1 can be justified for closed neural environments, in whichneurotransmitters cannot leave the synapse (no spill-over toother synapses). A2 is reasonable as long as only the so-called postsynaptic density [18] is considered, the part of thepostsynaptic membrane which contains most receptors. A3 canbe justified for low molecule concentrations. Finally, A4 is notobvious. For irreversible adsorption to a otherwise reflectivesurface covered by partially adsorbing disks, boundary homog-enization has been justified in [19] and refined with resultsfrom computer simulation in [20]. However, it is not clearif similar techniques can be applied for reversible reactions,too. There are two main issues to consider here. First, thesteady-state fluxes are substantially different; the net flux at areversibly adsorbing boundary is at steady-state, while forirreversible adsorption it is nonzero. Second, desorption altersthe spatial concentration profile of particles near the boundary;particles are more concentrated near receptors compared toirreversible adsorption [21]. The first issue can be resolved bycalibrating the effective adsorption coefficient to the homoge-nized surface such that it correctly reproduces the steady-stateflux to the patchy surface (see Section IV), while the secondissue can be resolved in a biologically plausible manner byassuming that particles may not re-adsorb immediately afterunbinding.With these assumptions, h ( t ) becomes independent of theparticle distribution in y and z , and, instead of C Ω ( x, y, z, t ) ,it is sufficient to consider the concentration of moleculesaggregated over y and z , C ( x, t ) , in µ m − , where C ( x, t ) = y max (cid:90) y min z max (cid:90) z min C Ω ( x, ν, ξ, t )d ξ d ν. (2) ig. 1. Model synapse. Neurotransmitters enclosed in vesicles at the presynaptic neuron are released into the synaptic cleft, propagate by Brownian motion,and activate receptors at the postsynaptic neuron. Binding to postsynaptic receptors is reversible. Furthermore, particles are re-uptaken (and recycled) at thepresynaptic neuron. The cuboid represents an abstraction from the synaptic cleft. Its orange and green surfaces represent the membrane of the pre- andpostsynaptic neuron, respectively. According to Fick’s second law of diffusion, (average)Brownian particle motion can be described by the followingpartial differential equation: ∂C∂t = D ∂ C∂x , (3)where D denotes the particle diffusion coefficient in µ m µ s − . B. Pre- and Postsynaptic Neurons
Particle re-uptake at the presynaptic neuron can be modeledas irreversible adsorption to the homogeneous left boundaryat x = x min (possibly after appropriate boundary homoge-nization) with re-uptake coefficient κ r in µ m µ s − . Accordingto the classical result in [22], adsorption of particles to apartially adsorbing boundary can be modeled with a radiatingboundary condition. To simplify notation, we set x min = 0 and x max = a . Then, the radiating boundary condition modelingpresynaptic re-uptake is given as D ∂C∂x = κ r C, at x = 0 . (4)For the right boundary, the radiating boundary needs tobe extended to incorporate particle desorption. We followa similar approach as described in [23]. Namely, particledesorption is modeled as a first-order process depending onthe intrinsic desorption rate κ d and on the amount of currentlyadsorbed particles, which, in turn, corresponds to h ( t ) . Thus,the boundary condition for reversible adsorption at the rightboundary can be formulated as D ∂C∂x = − adsorption (cid:122)(cid:125)(cid:124)(cid:123) κ a C + desorption (cid:122) (cid:125)(cid:124) (cid:123) κ d h ( t ) , at x = a (5) d h ( t )d t = κ a C ( a, t ) − κ d h ( t ) , (6)which simplifies directly to D ∂C∂x = − κ a C − κ d (cid:90) t D ∂C∂x d τ, at x = a. (7)Finally, instantaneous release of N particles at t = 0 ismodeled with the initial value C ( x,
0) =
N δ ( x − x ) , ≤ x ≤ a, (8) where δ ( x ) denotes the Dirac delta function. For the followingderivations, we set w.l.o.g. N = 1 . C. Channel Impulse Response
Once the solution to (3), (4), (7), (8) is found, h ( t ) can beobtained as h ( t ) = t (cid:90) − D ∂C ( x, τ ) ∂x (cid:12)(cid:12)(cid:12)(cid:12) x = a d τ. (9)III. A NALYTICAL C HANNEL M ODEL
A. Molecule Concentration
To find the C ( x, t ) that satisfies (3), (4), (7), (8), we use asimilar approach as [24, Ch. 14] and decompose C as C = U + W, (10)such that U fulfills (3) and (8), W fulfills (3) and equals at t , and U and W together fulfill (4) and (7). Next, weassume that the Laplace transform of C ( x, t ) with respect to(w.r.t.) t , ¯ C ( x, p ) = L{ C ( x, t ) } = (cid:82) ∞ C ( x, τ ) exp( − pτ )d τ ,exists. With this assumption and (10), C ( x, t ) can be found. Proposition 1:
Let κ r , κ d ≥ , κ a > . The unique solutionto (3), (4), (7), (8) is C ( x, t ) = ∞ (cid:88) n =1 Z n ( x ) Z n ( x ) e − Dα n t + κ d κ a + aκ d ( κ r = 0) , (11)where ( · ) denotes the indicator function, Z n ( x ) = [2 (cid:0) α n (cid:0) κ a − Dκ d (cid:1) + κ d + D α n (cid:1) ] D × ( Dα n cos ( xα n ) + κ r sin ( xα n )) , (12) D = a (cid:0) D α n + κ r (cid:1) (cid:0) α n (cid:0) κ a − Dκ d (cid:1) + κ d + D α n (cid:1) + Dα n ( Dκ d ( κ a − κ r ) + κ a κ r ( κ a + κ r ))+ κ d κ r ( κ a κ r + Dκ d ) + D α n ( κ a + κ r ) , (13) For κ a = 0 , h ( t ) = 0 , therefore this case is not considered here. nd the α n are defined as the positive roots of tan( aα ) = D ( κ a + κ r ) α − κ d κ r D α − ( κ a κ r + Dκ d ) α . (14) Proof:
Please refer to the Appendix.For irreversible adsorption to the right boundary (i.e., κ d =0 ), (11) reduces to [24, Ch. 14.3, eq. (4)], from which, inturn, follows [17, eq. (5)] after multiplying with the Green’sfunction for free diffusion in y and z . From now on, we assume κ d > .If there is no re-uptake, i.e., κ r = 0 , the steady-stateconcentration for large t is κ d κ a + aκ d , (15)because all exponential terms in (11) vanish. This is intuitive,as in the absence of re-uptake, particles are not ultimatelyremoved from the system. Also, the steady-state is a well-mixed state with constant concentration everywhere, indepen-dent of x . The net flux at the right boundary in this case iszero, meaning at each time step, the same number of particlesadsorb and desorb. Integrating over the size of the cleft, a , thefraction of solute particles in the steady-state is then given as aκ d κ a + aκ d , (16)while the fraction of permanently adsorbed particles is givenas − aκ d κ a + aκ d = κ a κ a + aκ d . (17)This result is a special case of [25, eq. (2.5)]. If aκ d (cid:29) κ a ,(16) approaches , i.e., almost all particles are solute, while,if κ a (cid:29) aκ d , the concentration of solute particles approaches and almost all particles are bound to receptors in the steady-state.In the presence of re-uptake (i.e. κ r > ), in contrast, thereis no constant term in (20), meaning that the concentrationof particles for t → ∞ approaches everywhere. Again, thisis intuitive, as in the presence of re-uptake, all particles areeventually re-uptaken. B. Channel Impulse Response
To obtain h ( t ) , we differentiate (11) with respect to x , andreach − ∂C ( x, t ) ∂x = ∞ (cid:88) n =1 Z (cid:48) n ( x ) Z n ( x ) e − Dα n t , (18)where Z (cid:48) n ( x ) = [2 (cid:0) α n (cid:0) κ a − Dκ d (cid:1) + κ d + D α n (cid:1) ] D × (cid:0) Dα n sin ( xα n ) − κ r α n cos ( xα n ) (cid:1) . (19)Evaluating (19) at x = a , assuming we may interchangeintegration and summation in (9), and computing the integralyields for the CIR h ( t ) = ∞ (cid:88) n =1 Z (cid:48) n ( a ) Z n ( x ) (cid:16) − e − Dα n t (cid:17) α n . (20) Due to the differentiation in (9), this term is independent ofthe constant term in (20) and thus valid for all κ r ≥ .We note that the sequence ( α n ) n is strictly monotonicallyincreasing. Therefore, for large t , the contributions of large- n terms are (almost) constant. In fact, the tail of h ( t ) can beproperly approximated with h ( t ) = − Z (cid:48) ( a ) Z ( x ) e − Dα t α , (21)see also Figs. 2–4. Such an approximation can be useful fordetector design and also a first step in quantifying ISI. How-ever, due to space constraints, further analytical investigationsare left for future work.IV. P ARTICLE - BASED S IMULATION
To verify the analytical expression derived in Section III,3-dimensional particle-based computer simulations were con-ducted. To this end, we adopted the simulator design from[26].
A. Simulator Design
Here, Brownian particle motion is simulated by updating theposition of each particle at each time step with a 3-dimensionaljointly independent Gaussian random vector [ X, Y, Z ] ∼ N ( × , σ I × ) , (22)where N ( µ , C ) denotes a multivariate Gaussian distributionwith mean vector µ and covariance matrix C , and × M and I M × M denote the × M all-zero vector and the M × M iden-tity matrix, respectively. Particle collisions are neglected. Thevariance σ is a function of the particle diffusion coefficient D and the simulation time step ∆ t , σ = 2 D ∆ t. (23)Accordingly, the root mean step length (rms) of a simulatedparticle is defined as s = √ D ∆ t. (24)For the simulation, the probability that a particle is re-uptaken after crossing the left boundary was computed from κ r using [26, eq. (21)]. Particles hitting a receptor at the postsy-naptic boundary were absorbed with probabilities computed asin [26, eqs. (37), (32)] from the intrinsic association coefficientof molecules to receptors, κ a , and the intrinsic desorption rateconstant, κ d . B. Boundary Homogenization
In order to compare simulation data and analytical results,boundary homogenization at the right boundary needs to beperformed. To avoid issues that arise from the non-uniformdistribution of desorbed particles, we choose ∆ t such that s islarger than the receptor radius, r . For fixed surface coverageat the postsynaptic neuron, increasing ∆ t has the effect thatdesorbing particles are more likely to see a representative partof the boundary before possible adsorbing to a receptor again.This is in fact equivalent to blocking the desorbed particlefor some time from re-adsorbing and thus fulfills one of theconditions of A4. Now, to compute the effective adsorption ABLE IS
IMULATION PARAMETERS FOR PARTICLE - BASED SIMULATION [27], [28].Parameter Default Value Description D . × − µ m µ s − Particle diffusion coefficient ∆ t Simulation time step N Number of released particles a
20 nm
Channel width in x direction y max − y min . µ m Channel width in y direction z max − z min . µ m Channel width in z direction κ (cid:48) r . Reduced re-uptake coefficient κ (cid:48) a . Reduced intrinsic adsorptioncoefficient κ (cid:48) d . Reduced desorption rate r . Receptor radius ρ . Receptor coverage at postsy-naptic neuron coefficient for the homogenized boundary, κ a , from κ a , weperform computer-assisted boundary homogenization similarto [20], with the difference, that we do not require analyticaland numerical results to produce the same average particlelife times , but instead demand matching steady-state concen-trations . To this end, we conduct particle-based simulationswithout re-uptake and let them run into steady-state. Next,we use our analytical result for the steady-state number ofadsorbed particles in the absence of re-uptake, (17), to fit κ a . In contrast to [20] and earlier results, it turns out that κ a for reversible reaction varies only moderately with r , if r is close to s , and mostly depends on the fraction of thepostsynaptic surface covered by receptors, ρ , and the intrinsicreceptor association coefficient κ a . We found that κ a = 0 . ρκ a (25)provides a good approximation for receptors of radius r =0 . , which is slightly less than s for ∆ t = 1 ns and D =6 . × − µ m µ s − .V. S IMULATION R ESULTS
To simplify the interpretation of the simulation resultspresented in this section, some of the simulation parametersare given as dimensionless quantities [26]. To this end, wedefine the reduced re-uptake coefficient κ (cid:48) r [26, eq. (12)] as κ (cid:48) r = κ r √ ∆ t √ D , (26)the reduced intrinsic adsorption coefficient κ (cid:48) a [26, eq. (12)]as κ (cid:48) a = κ a √ ∆ t √ D , (27)and, finally, the reduced desorption rate κ (cid:48) d [26, eq. (13)] as κ (cid:48) d = κ d ∆ t. (28)For the diffusion coefficient and the width of the synapticcleft, we used values from the literature. The other parameterswere varied to ensure the validity of our model for a widerange of sensible parameter values. A complete listing of thedefault parameters can be found in Table I.As particles diffuse independently, the CIR for N > is obtained by multiplying (20) by N . For the numerical Fig. 2. Results for particle-based simulation (gray curves) and analyticalsolutions (solid lines) for different receptor densities at the postsynapticneuron. The first-term approximation (21) is shown with a dashed line. evaluation of (20), the sum was truncated after the first terms. After proper calibration of the homogenized adsorptioncoefficient, κ a , as described in Section IV, the agreementbetween analytical solution and simulation results was ingeneral excellent for all parameter regimes that were tested.First, we investigate the impact of the receptor coverageat the postsynaptic neuron, ρ , on the CIR. In Fig. 2, itcan be seen that, on the one hand, increasing the receptorcoverage increases the peak value, but, on the other hand,the CIR also takes more time to decay. This is expected asby increasing the receptor coverage, firstly, more particlesadsorb to the postsynaptic membrane at their first arrival,and, secondly, desorbed particles are more likely to rebindimmediately after desorption. Also, Fig. 2 shows that the CIRcan be approximated with high accuracy after its peak usingthe first term of the sum in (20).Next, we investigate the effect of particle re-uptake on theCIR. From Fig. 3, it can be observed that increasing the re-uptake rate considerably shortens the CIR. At the same time,however, the peak value of the received signal is reduced.Again, this trade-off is expected as, for larger κ (cid:48) r , particlesare more likely to be re-uptaken before they get the chance toeven reach the postsynaptic side.Finally, the impact of the channel width on the CIR isinvestigated in Fig. 4. Increasing the channel width leads toa more dispersive CIR and a more severe attenuation of thesignal. Interestingly, in contrast to the other parameters weinvestigated in Figs. 2 and 3, here, decreasing the cleft widthleads to a shorter and stronger signal. Thus, in terms of ISImitigation and peak value, a small cleft width is beneficial.VI. C ONCLUSIONS
In this paper, a new analytical model for the synapticchannel has been proposed. An analytical time domain solutionfor the CIR has been derived and validated with particle-based simulation. The supposed effect of presynaptic moleculere-uptake, namely a significant shortening of the CIR, hasbeen observed in the model system. Furthermore, it has beendemonstrated that molecule re-uptake also leads to a lowerpeak value. It was shown that higher receptor density has a ig. 3. Results for particle-based simulation (gray curves) and analyticalsolutions (solid lines) are shown for different re-uptake rates κ r togetherwith the first-term approximation (21) (dashed lines).Fig. 4. Results for particle-based simulation (gray curves) and analyticalsolutions (solid lines) are shown for different cleft widths together with thefirst-term approximation (21) (dashed lines). positive effect on the peak value, while, at the same time,also the undesired signal part (tail) is enhanced. In addition,the impact of the width of the synaptic cleft on the CIR wasinvestigated and it was shown that a smaller distance betweentransmitter and receiver is beneficial. Finally, the first term ofthe infinite sum in the proposed CIR expression was shown toprovide an easy-to-compute approximation for the tail of theCIR.We believe that by incorporating many important bio-physical features of the synaptic communication channel, theproposed model is useful for the design of future syntheticneural communication systems.A PPENDIX
A. Sketch of Derivation of CIR
Eq. (4) reads in the Laplace domain
D ∂ ¯ C∂x = κ r ¯ C, for x = 0 , (29)while (7) transforms to D ∂ ¯ C∂x = − κ a ¯ C − κ d D p ∂ ¯ C∂x , for x = a. (30) Here, p is the Laplace variable and we define q = (cid:112) pD .Let ¯ U and ¯ W denote the Laplace transforms of U and W ,as given in (10), respectively.Following [24, Ch. 14.3], ¯ U and ¯ W can be found to be ¯ U = 12 Dq e − q | x − x (cid:48) | , (31)and ¯ W = A sinh( qx ) + B cosh( qx ) , (32)where constants A and B are to be chosen such that (29) and(30) are fulfilled.Substituting ¯ C = ¯ U + ¯ W in (29) and (30) and solving for A and B , we obtain ¯ C as ¯ C ( x, p ) = e − q ( a + | x − x | ) (cid:16) Dqe q ( a + | x − x |− x ) (cid:0) cosh( q ( a − x )) (cid:0) κ d + Dq (cid:1) + qκ a sinh( q ( a − x ))) − κ d κ r cosh( q ( a − x )) e q ( a + | x − x |− x ) − Dq κ r cosh( q ( a − x )) e q ( a + | x − x |− x ) − Dq κ a cosh( qx ) e q ( | x − x | + x ) − qκ a κ r sinh( q ( a − x )) e q ( a + | x − x |− x ) − qκ a κ r sinh( qx ) e q ( | x − x | + x ) + Dqe aq κ d sinh( aq )+ e aq κ d κ r cosh( aq ) + D q e aq sinh( aq )+ Dq κ a e aq cosh( aq ) + Dq e aq κ r cosh( aq )+ qκ a e aq κ r sinh( aq )+ e q ( | x − x | + x ) (cid:0) κ d + Dq (cid:1) ( Dq cosh( qx ) + κ r sinh( qx )) (cid:17) / (cid:2) Dq (cid:0) q sinh( aq ) (cid:0) κ a κ r + D (cid:0) κ d + Dq (cid:1)(cid:1) + cosh( aq ) (cid:0) Dq κ a + κ r (cid:0) κ d + Dq (cid:1)(cid:1)(cid:1)(cid:3) . (33) B. Time Domain Solution
The corresponding solution in the time domain is now givenby the inverse Laplace transform: C ( x, t ) = 12 πj (cid:90) γ + j ∞ γ − j ∞ e pt ¯ C ( x, p ) dp, (34)where j denotes the imaginary unit and γ ∈ R + needs to bechosen large enough such that all singularities of ¯ C are leftof the line along which the integral is computed.Similarly to [24], instead of integrating along the infiniteline, we complete it with a semicircle to a closed contourwhich contains the origin and all singularities of ¯ C . Then, weuse the residue theorem to replace this contour integral witha sum over the residues { σ } of ¯ C , C ( x, t ) = (cid:88) { σ } Res( e pt ¯ C, σ ) . (35)First, we note that ¯ C can be written as quotient of twofunctions, f and g , which are given as the numerator anddenominator of (33), respectively. The residues of e pt ¯ C at allsimple poles σ can then be computed by l’Hˆopital’s rule asRes (cid:0) e pt ¯ C, σ (cid:1) = e σt f ( x, σ ) g (cid:48) ( x, σ ) . (36)For higher order poles, the residues can be found by expanding e pt ¯ C as Laurent series.et us now look at the denominator of ¯ C , Dq (cid:0) q sinh( aq ) (cid:0) κ a κ r + D (cid:0) κ d + Dq (cid:1)(cid:1) + cosh( aq ) (cid:0) Dq κ a + κ r (cid:0) κ d + Dq (cid:1)(cid:1)(cid:1) . (37)The parameters D and a are always positive.Now, let us first assume that κ r , κ a , and κ d are also positive.In this case, (37) has one simple pole at q = 0 and nonzerosimple poles at the roots of sinh( aq )cosh( aq ) = Dq κ a + κ r (cid:0) κ d + Dq (cid:1) q ( κ a κ r + D ( κ d + Dq )) . (38)Set q = jα , then the non-zero singularities of ¯ C are given asthe roots ± α n , n ∈ N , of tan( aα ) = D ( κ a + κ r ) α − κ d κ r D α − ( κ a κ r + Dκ d ) α . (39)Res (cid:0) e pt ¯ C, (cid:1) can be computed to beRes (cid:0) e pt ¯ C, (cid:1) = 02 Dκ r κ d = 0 . (40)By repeatedly exploiting (39) and finally plugging in p = − Dα n , the residues of e pt ¯ C at the simple non-zero poles α n can be computed to beRes (cid:0) e pt ¯ C, α n (cid:1) = 2 (cid:0) α n (cid:0) κ a − Dκ d (cid:1) + κ d + D α n (cid:1) D× ( Dα n cos ( xα n ) + κ r sin ( xα n )) × ( Dα n cos ( α n x (cid:48) ) + κ r sin ( α n x (cid:48) )) , (41)where D is defined in (13). Factorizing (41) and summing overall α n yields (11).Now, if κ a > and κ d > , but κ r = 0 , (37) has a doubleroot at q = 0 and Res (cid:0) e pt ¯ C, (cid:1) is most easily computed fromthe Laurent series expansion of e pt ¯ C . Namely, the residue of e pt ¯ C coincides with the coefficient a − of its Laurent seriesexpansion at p = 0 .Expanding (37) in q yields g ( x,
0) = q (2 aDκ d + 2 Dκ a ) + O (cid:0) q (cid:1) . (42)Because p = Dq , the coefficient of p − in the expansion of e pt ¯ C corresponds to the coefficient of ( Dq ) − and it is clearthat only coefficients of constant terms from the expansionof the numerator in q play a role. Thus, expanding the terms κ d cosh( q ( a − x )) and κ d cosh( qx ) and dividing by (42) yieldsRes (cid:0) e pt ¯ C, (cid:1) = κ d κ a + aκ d . (43)For κ d = 0 , C ( x, t ) can be obtained in a similar fashion using(36). This completes the proof.R EFERENCES[1] T. Nakano, A. W. Eckford, and T. Haraguchi,
Molecular Communica-tion . Cambridge University Press, 2013.[2] 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.[3] R. A. Freitas,
Nanomedicine, Volume I: Basic Capabilities . LandesBioscience Georgetown, TX, 1999, vol. 1.[4] M. Veleti´c and I. Balasingham, “Synaptic communication engineeringfor future cognitive brain–machine interfaces,”
Proc. IEEE , vol. 107,no. 7, pp. 1425–1441, Jul. 2019. [5] D. W. K. Ng, E. S. Lo, and R. Schober, “Wireless information andpower transfer: Energy efficiency optimization in OFDMA systems,”
IEEE Trans. Wireless Commun. , vol. 12, no. 12, pp. 6352–6370, 2013.[6] Y. Deng, W. Guo, A. Noel, A. Nallanathan, and M. Elkashlan, “Enablingenergy efficient molecular communication via molecule energy transfer,”
IEEE Commun. Lett. , vol. 21, no. 2, pp. 254–257, Feb. 2017.[7] W. Guo, Y. Deng, H. B. Yilmaz, N. Farsad, M. Elkashlan, A. Eckford,A. Nallanathan, and C. Chae, “SMIET: Simultaneous molecular infor-mation and energy transfer,”
IEEE Wireless Commun. , vol. 25, no. 1,pp. 106–113, Feb. 2018.[8] H. Arjmandi, M. Movahednasab, A. Gohari, M. Mirmohseni, M. Nasiri-Kenari, and F. Fekri, “ISI-avoiding modulation for diffusion-basedmolecular communication,”
IEEE Trans. Mol. Biol. Multi-Scale Com-mun. , vol. 3, no. 1, pp. 48–59, Mar. 2017.[9] M. S. Leeson and M. D. Higgins, “Forward error correction for molecu-lar communications,”
Nano Commun. Networks , vol. 3, no. 3, pp. 161–167, 2012.[10] P. Shih, C. Lee, P. Yeh, and K. Chen, “Channel codes for reliability en-hancement in molecular communication,”
IEEE J. Sel. Areas Commun. ,vol. 31, no. 12, pp. 857–867, Dec. 2013.[11] B. Tepekule, A. E. Pusane, H. B. Yilmaz, C. Chae, and T. Tugcu, “ISImitigation techniques in molecular communication,”
IEEE Trans. Mol.Biol. Multi-Scale Commun. , vol. 1, no. 2, pp. 202–216, Jun. 2015.[12] A. Noel, K. C. Cheung, and R. Schober, “Improving receiver per-formance of diffusive molecular communication with enzymes,”
IEEETrans. Nanobiosci. , vol. 13, no. 1, pp. 31–43, Mar. 2014.[13] A. C. Heren, H. B. Yilmaz, C. Chae, and T. Tugcu, “Effect of degrada-tion in molecular communication: Impairment or enhancement?”
IEEETrans. Mol. Biol. Multi-Scale Commun. , vol. 1, no. 2, pp. 217–229, Jun.2015.[14] A. S. Kristensen, J. Andersen, T. N. Jørgensen, L. Sørensen, J. Eriksen,C. J. Loland, K. Strømgaard, and U. Gether, “SLC6 neurotransmit-ter transporters: Structure, function, and regulation,”
Pharmacol. Rev. ,vol. 63, no. 3, pp. 585–640, 2011.[15] E. Balevi and O. B. Akan, “A physical channel model for nanoscaleneuro-spike communications,”
IEEE Trans. Commun. , vol. 61, no. 3,pp. 1178–1187, Mar. 2013.[16] M. Veleti´c, F. Mesiti, P. A. Floor, and I. Balasingham, “Communicationtheory aspects of synaptic transmission,” in
Proc. IEEE Intern. Conf.Commun. , 2015, pp. 1116–1121.[17] T. Khan, B. A. Bilgin, and O. B. Akan, “Diffusion-based model forsynaptic molecular communication channel,”
IEEE Trans. Nanobiosci. ,vol. 16, no. 4, pp. 299–308, Jun. 2017.[18] M. Sheng and C. C. Hoogenraad, “The postsynaptic architecture ofexcitatory synapses: A more quantitative view,”
Annu. Rev. Biochem. ,vol. 76, no. 1, pp. 823–847, 2007, pMID: 17243894.[19] R. Zwanzig and A. Szabo, “Time dependent rate of diffusion-influencedligand binding to receptors on cell surfaces,”
Biophys. J. , vol. 60, no. 3,pp. 671–678, 1991.[20] A. M. Berezhkovskii, Y. A. Makhnovskii, M. I. Monine, V. Y. Zitserman,and S. Y. Shvartsman, “Boundary homogenization for trapping by patchysurfaces,”
J. Chem. Phys. , vol. 121, no. 22, pp. 11 390–11 394, 2004.[21] A. Szabo, “Theoretical approaches to reversible diffusion-influencedreactions: Monomer–excimer kinetics,”
J. Chem. Phys. , vol. 95, no. 4,pp. 2481–2490, 1991.[22] F. C. Collins and G. E. Kimball, “Diffusion-controlled reaction rates,”
J. Colloid Sci. , vol. 4, no. 4, pp. 425–437, 1949.[23] A. Ahmadzadeh, H. Arjmandi, A. Burkovski, and R. Schober, “Com-prehensive reactive receiver modeling for diffusive molecular commu-nication systems: Reversible binding, molecule degradation, and finitenumber of receptors,”
IEEE Trans. Nanobiosci. , vol. 15, no. 7, pp. 713–727, Oct. 2016.[24] H. Carslaw and J. Jaeger,
Conduction of Heat in Solids , ser. Oxfordscience publications. Clarendon Press, 1986.[25] A. M. Berezhkovskii and A. Szabo, “Effect of ligand diffusion onoccupancy fluctuations of cell-surface receptors,”
J. Chem. Phys. , vol.139, no. 12, p. 121910, 2013.[26] S. S. Andrews, “Accurate particle-based simulation of adsorption, des-orption and partial transmission,”
Phys. Biol. , vol. 6, no. 4, p. 046015,Nov. 2009.[27] M. Rice, G. Gerhardt, P. Hierl, G. Nagy, and R. Adams, “Diffusion coef-ficients of neurotransmitters and their metabolites in brain extracellularfluid space,”
Neuroscience , vol. 15, no. 3, pp. 891–902, 1985.[28] B. Alberts, D. Bray, K. Hopkin, A. Johnson, J. Lewis, M. Raff,K. Roberts, and P. Walter,