We set up a coordinate system such that the *x* and *y* axes are parallel to the pre- and post-synaptic membranes. The synaptic junction has depth *L*_{
z
} . The junction is periodic in the *x* and *y* directions, with periodicities of *L*_{
x
} and *L*_{
y
} , respectively. In each "unit cell", a synaptic vesicle bursts at time *t* = 0, releasing the neurotransmitters into the cleft. We model the release of the neurotransmitters as a transient flux, *u*(*t*), that is confined to a circular opening with radius *R*. We place the synaptic vesicle at the center of the pre-synaptic face of the unit cell. After diffusing to the post-synaptic membrane, neurotransmitters are absorbed by a circular disk in each unit cell, with radius *a*. We place this "sink" also at the center of the post-synaptic face of the unit cell. The exact shapes and locations of the pre-synaptic opening and the post-synaptic sink are not essential for the analytical solution of our model. The quantity of interest is the total flux, *J*(*t*), through the post-synaptic face of each unit cell.

We place the origin of the Cartesian coordinate system at the center of the pre-synaptic face, with the *z* axis pointing toward the post-synaptic face. The concentration of neurotransmitters at position **r** and at time *t* is *C*(**r**, *t*). Within the synaptic junction, it is governed by the diffusion equation,

\frac{\partial C(\mathbf{r},t)}{\partial t}=D{\nabla}^{2}C(\mathbf{r},t)

(1)

where *D* is the diffusion constant. The initial condition is

The boundary condition at *z* = 0 is

D\frac{\partial C(\mathbf{r},t)}{\partial z}=u(t)\phantom{\rule{0.5em}{0ex}}\text{when}\phantom{\rule{0.5em}{0ex}}\rho <R

(3a)

D\frac{\partial C(\mathbf{r},t)}{\partial z}=0\phantom{\rule{0.5em}{0ex}}\text{when}\phantom{\rule{0.5em}{0ex}}\rho >R

(3b)

where *ρ* = (*x*^{2} + *y*^{2})^{1/2} is the distance to the *z* axis. The boundary condition at *z* = *L*_{
z
} is

C(\mathbf{r},t)=0\phantom{\rule{0.5em}{0ex}}\text{when}\phantom{\rule{0.5em}{0ex}}\rho <a

(4a)

D\frac{\partial C(\mathbf{r},t)}{\partial z}=0\phantom{\rule{0.5em}{0ex}}\text{when}\phantom{\rule{0.5em}{0ex}}\rho >a

(4b)

We solve the problem in Lapace space. For a function *f*(*t*) of *t*, we denote the Laplace transform as \widehat{f}(s)\equiv {\displaystyle \underset{0}{\overset{\infty}{\int}}dt}{e}^{-st}f(t). The Laplace transform of Eq. (1), using the initial condition of Eq. (2), is

s\widehat{C}(\mathbf{r},s)=D{\nabla}^{2}\widehat{C}(\mathbf{r},s)

(5)

The solution appropriate for the periodic boundary conditions in the *x* and *y* directions has the form

\widehat{C}(\mathbf{r},s)={\displaystyle \sum _{l,m=0}^{\infty}\left[{\alpha}_{lm}(s){e}^{{\gamma}_{lm}(s)z}+{\beta}_{lm}(s){e}^{-{\gamma}_{lm}(s)z}\right]\mathrm{cos}(2l\pi x/{L}_{x})\mathrm{cos}(2m\pi y/{L}_{y})}

(6)

where

{\gamma}_{lm}(s)={\left[s/D+{(2l\pi /{L}_{x})}^{2}+{(2m\pi y/{L}_{y})}^{2}\right]}^{1/2}

(7)

The coefficients *α*_{
lm
} (*s*) and *β*_{
lm
} (*s*) are to be determined from boundary conditions.

To make use of the boundary condition of Eq. (3), we need the 2-dimensional cosine transform of a function in *x* and *y* that has value 1 when *ρ* < *R* and value 0 when *ρ* > *R*. The coefficient of the cos(2*lπ*/*L*_{
x
} )cos(2*mπ*/*L*_{
y
} ) term is

{q}_{lm}=\frac{4{\epsilon}_{l}{\epsilon}_{m}}{{L}_{x}{L}_{y}}{\displaystyle \underset{\rho <R}{\iint}dxdy}\mathrm{cos}(2l\pi x/{L}_{x})\mathrm{cos}(2m\pi y/{L}_{y})

(8)

where *ε*_{0} = 1/2 and *ε*_{
l
} = 1 for *l* > 0. For later reference, we note that the function having value 1 when *ρ* < *a* and value 0 when *ρ* > *a* has the following coefficients in its 2-dimensional cosine transform:

{p}_{lm}=\frac{4{\epsilon}_{l}{\epsilon}_{m}}{{L}_{x}{L}_{y}}{\displaystyle \underset{\rho <a}{\iint}dxdy}\mathrm{cos}(2l\pi x/{L}_{x})\mathrm{cos}(2m\pi y/{L}_{y})

(9)

Applying the boundary condition of Eq. (3), we find

D{\gamma}_{lm}(s)\left[{\alpha}_{lm}(s)-{\beta}_{lm}(s)\right]=\stackrel{\u02c6}{u}(s){q}_{lm}

(10)

For the boundary condition of Eq. (4) at *z* = *L*_{
z
} , we use the constant-flux approximation [9]:

D\frac{\partial \widehat{C}(\mathbf{r},s)}{\partial z}=Q(s)\phantom{\rule{0.5em}{0ex}}\text{when}\phantom{\rule{0.5em}{0ex}}\rho <a

(11a)

D\frac{\partial \widehat{C}(\mathbf{r},s)}{\partial z}=0\phantom{\rule{0.5em}{0ex}}\text{when}\phantom{\rule{0.5em}{0ex}}\rho >a

(11b)

where the quantity *Q*(*s*) is to be determined by the condition

<\widehat{C}(\mathbf{r},s)>\equiv {(\pi {a}^{2})}^{-1}{\displaystyle \underset{\rho <a}{\iint}dxdy}\widehat{C}(\mathbf{r},s)=0

(12)

Equation (11) leads to

D{\gamma}_{lm}(s)\left[{\alpha}_{lm}(s){e}^{{\gamma}_{lm}(s){L}_{z}}-{\beta}_{lm}(s){e}^{-{\gamma}_{lm}(s){L}_{z}}\right]=Q(s){p}_{lm}

(13)

With Eqs. (10) and (13), we solve for coefficients *α*_{
lm
} (*s*) and *β*_{
lm
} (*s*). The results are

{\alpha}_{lm}{e}^{{\gamma}_{lm}(s){L}_{z}}=\frac{1}{D{\gamma}_{lm}(s)}\frac{-{q}_{lm}\stackrel{\u02c6}{u}(s)+{p}_{lm}Q(s){e}^{{\gamma}_{lm}(s){L}_{z}}}{{e}^{{\gamma}_{lm}(s){L}_{z}}-{e}^{-{\gamma}_{lm}(s){L}_{z}}}

(14a)

{\beta}_{lm}{e}^{-{\gamma}_{lm}(s){L}_{z}}=\frac{1}{D{\gamma}_{lm}(s)}\frac{-{q}_{lm}\stackrel{\u02c6}{u}(s)+{p}_{lm}Q(s){e}^{-{\gamma}_{lm}(s){L}_{z}}}{{e}^{{\gamma}_{lm}(s){L}_{z}}-{e}^{-{\gamma}_{lm}(s){L}_{z}}}

(14b)

Inserting Eq. (6) in Eq. (12) and using Eqs. (14) and (9), we find *Q*(*s*) as

Q(s)=\stackrel{\u02c6}{u}(s)\frac{{\displaystyle \sum _{l,m=0}^{\infty}\frac{{p}_{lm}{q}_{lm}}{{\epsilon}_{l}{\epsilon}_{m}{\gamma}_{lm}(s)\mathrm{sinh}[{\gamma}_{lm}(s){L}_{z}]}}}{{\displaystyle \sum _{l,m=0}^{\infty}\frac{{p}_{lm}{}^{2}\mathrm{cosh}[{\gamma}_{lm}(s){L}_{z}]}{{\epsilon}_{l}{\epsilon}_{m}{\gamma}_{lm}(s)\mathrm{sinh}[{\gamma}_{lm}(s){L}_{z}]}}}

(15)

Finally the total flux through the sink on the post-synaptic face is

\widehat{J}(s)=\pi {a}^{2}Q(s)

(16)

The flux accumulated over all times,

I={\displaystyle \underset{0}{\overset{\infty}{\int}}dtJ(t)}

(17a)

is of interest. In our model,

I=\widehat{J}(0)=\pi {R}^{2}\widehat{u}(0)

(17b)

which is independent of *a*. Equation (17b) is simply a consequence of ligand conservation, i.e., the total number of ligands released from the synaptic vesicle is the same as the total number of ligands absorbed by the receptors.

The analytical solution derived above can be implemented on any function *u*(*t*) modeling neurotransmitter release from the synaptic vesicle. We focused on an exponentially decaying function:

u(t)={e}^{-t/{t}_{0}}

(18a)

Its Laplace transform is

\widehat{u}(s)=\frac{1}{s+1/{t}_{0}}

(18b)

Our model then contains two parameters related to time: *D* and *t*_{0}. For two sets of these parameters, e.g., *D*_{1} and *t*_{01} in one and *D*_{2} and *t*_{02} in the other, it can be shown that the corresponding response functions satisfy

{J}_{1}({D}_{2}t)={J}_{2}({D}_{1}t)\phantom{\rule{0.5em}{0ex}}\text{when}\phantom{\rule{0.5em}{0ex}}{D}_{1}{t}_{01}={D}_{2}{t}_{02}

(19)

We now briefly describe the details of our implementation of the analytical solution. To calculate the *q*_{
lm
} and *p*_{
lm
} coefficients of Eqs. (8) and (9), we first carried out the integration over *y* analytically. The remaining integration over *x* was done numerically using the Gauss-Legendre quadrature with two points. The summations over *l* and *m* in Eq. (15) were truncated at *l* = *m* = 40. The Laplace transform of \widehat{J}(s) was inverted by the Stehfest algorithm [10]. A Fortran90 code for the implementation is available upon request.