Back Action on Neurotransmitters by Receptor Binding Reveals an Optimal Receptor Density Profile

In this study, we illustrate the effect of such back action on structural optimization by focusing on a particular example of postsynaptic receptors. Transmission of information through chemical synapses consists of three parts: release of chemical messengers (neurotransmitters) from one neuron into the extracellular space (synaptic cleft), neurotransmitter diffusion in the synaptic cleft, and neurotransmitter binding to receptors located on the receiving (postsynaptic) neuron. The process of neurotransmitter binding and unbinding (to receptors) is usually assumed not to modify the concentration of free neurotransmitter in the synaptic cleft, although this assumption is only valid when the number of neurotransmitter molecules is very high relative to the number of receptors and the diffusion rate is fast relative to the binding kinetics, which is not necessarily the case. How the dynamics of diffusing ligands are altered in the presence of interactions (such as binding/unbinding to ligand binding sites or consumption via oxidation) has been extensively studied both experimentally and theoretically. Very early examples involve the study of oxygen diffusion through tissue and its consumption via oxidation [1,2]. Early experiments [3-5] involving neurotransmitter diffusion indicated that the binding of neurotransmitters to receptors can effectively slow the diffusion of neurotransmitters, a phenomenon often referred to as buffered diffusion, especially when binding kinetics are much faster than diffusion rates, and additional experiments [6] have confirmed other predictions of buffered diffusion models. There have also been numerous theoretical works, involving deterministic studies using partial differential equations [7–11] and stochastic studies using Monte Carlo simulations [12,13], which have revealed the importance of incorporating receptor-neurotransmitter interactions in evaluating free neurotransmitter concentration. However, most of these studies have focused on the changes in the temporal characteristics of ligand-receptor interactions due to buffered diffusion rather than on its resulting functional implications which we investigate in this study.


Introduction
Coupled rate equations are commonly used to describe biophysical time evolution processes on a high course-grained level. The rate constants are typically treated as fitting parameters, which contain information on the microscopic details underlying the dynamical processes. Typically, the rate equations form a differential linear system of equations given by ( ) ( ) ( ) = A   d y t t y t dt (1) Where external sources are treated as being unperturbed by the evolution of the variables ( )  y t and are absorbed into the matrix A(t). In this work we consider a generalization of this approach, which includes the effect of the back action of the variables ( )  y t on the source terms. The resulting set of coupled differential equations is non-linear, and it encodes temporal feedback of the system on its environment.
In this study, we illustrate the effect of such back action on structural optimization by focusing on a particular example of postsynaptic receptors. Transmission of information through chemical synapses consists of three parts: release of chemical messengers (neurotransmitters) from one neuron into the extracellular space (synaptic cleft), neurotransmitter diffusion in the synaptic cleft, and neurotransmitter binding to receptors located on the receiving (postsynaptic) neuron. The process of neurotransmitter binding and unbinding (to receptors) is usually assumed not to modify the concentration of free neurotransmitter in the synaptic cleft, although this assumption is only valid when the number of neurotransmitter molecules is very high relative to the number of receptors and the diffusion rate is fast relative to the binding kinetics, which is not necessarily the case. How the dynamics of diffusing ligands are altered in the presence of interactions (such as binding/unbinding to ligand binding sites or consumption via oxidation) has been extensively studied both experimentally and theoretically. Very early examples involve the study of oxygen diffusion through tissue and its consumption via oxidation [1,2]. Early experiments [3][4][5] involving neurotransmitter diffusion indicated that the binding of neurotransmitters to receptors can effectively slow the diffusion of neurotransmitters, a phenomenon often referred to as buffered diffusion, especially when binding kinetics are much faster than diffusion rates, and additional experiments [6] have confirmed other pre-dictions of buffered diffusion models. There have also been numerous theoretical works, involving deterministic studies using partial differential equations [7][8][9][10][11] and stochastic studies using Monte Carlo simulations [12,13], which have revealed the importance of incorporating receptor-neurotransmitter interactions in evaluating free neurotransmitter concentration. However, most of these studies have focused on the changes in the temporal characteristics of ligand-receptor interactions due to buffered diffusion rather than on its resulting functional implications which we investigate in this study.
Incorporating receptor-neurotransmitter interactions is expected to have a significant effect when the number of receptors can be much larger than the number of ligand molecules and the kinetics of binding/ unbinding can be very fast relative to the free diffusion rate [5]. On the other hand, when the number of neurotransmitter molecules far exceeds the number of receptors and/or when binding is slow relative to diffusion rates, the effects of binding/unbinding are expected to be minimal. In these cases, it is generally assumed that the concentration of free neurotransmitter is not significantly modified by the binding of the neurotransmitter to the receptors. For instance, if the number of postsynaptic glutamate receptors is small relative to the number of glutamate molecules being released by a presynaptic event, it is logical to expect that the binding/unbinding of glutamate to its receptors will minimally affect free glutamate concentration [14,15] where binding/ unbinding is neglected in the modelling). This assumption is further supported by experimental results indicating that a single vesicle releases approximately 3000 glutamate molecules, and that the number of postsynaptic glutamate receptors is on the order of 100 receptors [16]. It would therefore seem that the number of receptors is sufficiently small relative to the number of glutamate molecules to neglect the effect of glutamate binding to its receptor on free glutamate concentration. However, since the association rate constant of glutamate to AMPA receptors is of the same order of magnitude as the diffusion rate, glutamate binding to receptors could significantly decrease the number of glutamate molecules as the distance from the release site increases, suggesting that it is not obvious that this assumption is valid.
The goal of this study is to apply a rate equation approach to determine whether neurotransmitter/receptor interactions would be different with a free diffusion model compared to a receptor-ligand diffusion model. We will investigate the temporal characteristics of ligand-receptor interactions due to buffered diffusion as well as its resulting quantitative functional implications on global synaptic AMPA receptor function. We consider a classical receptor kinetic model for the AMPA receptors, which has been described in numerous publications [17][18][19] using a glutamate diffusion profile unaffected by binding (free dif-fusion case). We then determine the dynamics of the receptor-limited diffusion model, this time considering the impact of glutamate binding/unbinding on free glutamate concentration profile in the cleft. We discover that there is a clear difference between the free diffusion and the receptor-limited diffusion model, thereby underscoring that the free diffusion case, which is most commonly used, may lead to inaccuracies in the calculation of transmitter/ receptor interactions. Moreover, when we take into consideration actual dimensions of the postsynaptic densities, our results indicate the existence of an optimal receptor density for a given radial distance from the release site, a feature that cannot be accounted for with the free diffusion model. Interestingly, this optimal density is on the order of the experimentally determined value. These results suggest that the dimensions of postsynaptic densities and the receptor density are set by some basic physical properties of ligand-receptor interactions.

Glutamate-AMPA receptor model
An illustration of the glutamate/receptor system we are considering is shown in Figure 1. Glutamate is released from the presynaptic terminal in the synaptic cleft in the form of a disk of radius r Glu and height h, which corresponds respectively to the size of vesicle, and the distance between the pre-and post-synaptic elements (the transverse length of the cleft) [20,21]. Diffusion of glutamate takes place over a disk of receptors of radius r PSD . Glutamate molecules are considered to be small enough to be described in the continuum limit, and we assume they are uniformly distributed throughout the cylinder. We assume that the height h is sufficiently small for the dynamics along the height direction to be negligible, such that any changes in glutamate concentration occur instantaneously along this direction. This has been shown to be a valid approximation for typical cleft heights [22]. In turn, the only relevant diffusion occurs in the plane transverse to the height.
The dynamics of AMPA receptors are governed by the kinetic schema for the 16-state AMPA receptor model [17], depicted in Figure 2. The model describes the opening of ion channels due to the binding of Glutamate to the receptor. The states in the schema are labeled by R,D,E,O, where D refers to being singly desensitized, E to doubly desensitized, and O to the opening of the receptor-associated ion channel. Furthermore, states are followed by a number 1,2,3,4 referring to how many glutamate molecules are bound. The states take values from 0 to 1 corresponding to the fraction of the total number of receptors in that state. Moving to the right along the diagram requires the binding of glutamate (measured in mM). Although a simpler model could be used to study the phenomena of interest, we selected this model because results generated by this model provide a good fit with numerous experimental data under a variety of experimental conditions. Furthermore, as will be discussed below, the existence of multiple binding sites for each molecule of receptors plays a significant role in the difference between the 2 cases. The parameters for the model given in Table 1 were adapted from the original values used in Robert and Howe [17] to additionally fit the experimental results from Kessler et al. [23]. This optimization was presented in Bouteiller et al. [19]. We write the coupled rate equations associated with the kinetic schemes as follows:
where x is the spatial coordinate, t is time, and for conciseness the states in the schema are labelled by y i , with i=0,…, 15 and y 16 corresponds to glutamate concentration. We explain the form of the rate equation for glutamate and the function F (y(x, t)) below.

Rate equation for glutamate concentration
The number of glutamate molecules in the system is given by: Where ( ) j i y is the fraction of receptors in the state i with J glutamate molecules bound to receptors cleft N is the number of molecules in the cleft (which are not bound to a receptor) and N is the number of receptors. Since the dynamics we are considering do not create or destroy glutamate molecules, taking a time derivative of equation (3) gives zero on the left-hand side, which then gives us an equation for the time evolution of the molecules in the cleft: The derivative acting on the y i 's can be replaced with the coupled rate equations in equation (2). We can now convert this equation to describe the concentration of glutamate at a given location x as: where ρ rec is the density of receptors (in mM) associated with the spatial distribution of receptors. For notational conciseness, we define: For simplicity we assume a "smearing" of the receptors such that they can be described by a continuous distribution. For example, a uniform distribution of receptors in a disk of radius r PSD , ρ rec would be given by: Where h D is the Debye length (µm) Here, the Debye length corresponds to an interaction distance, representing the effective length along the height of the cylinder of Figure 1 from which the receptors can bind/unbind glutamate. To include diffusion, we simply include the diffusion term: Where D Glu is the (free) diffusion constant for glutamate in the cleft

Numerical methods
The system of partial differential equations (n+1 in total, where n=16) we have to solve is given by: Since the system and initial conditions are rotationally symmetric about the origin (this means we have Neumann boundary conditions for glutamate concentration at the origin), the angular dependence in the diffusion term can be ignored: Where r denotes the radial coordinate. We discretize r with a grid of m elements with a maximum radius of 1 μm, and then we consider the time evolution of the m×(n+1) dependent variables. Since we are discretizing in the radial direction, our grid is formed from annular slices of width where we take m=250 for our simulations. At the maximum radius of 1 μm, we impose Dirichlet boundary conditions (an absorbing boundary) on the glutamate concentration to ensure that any glutamate molecule that has diffused to this point is irrecoverable. To solve this problem, we use the numerical method TR-BDF2 with an adaptive step size [24,25].
We focus on the case where receptors have a uniform distribution within a disk of fixed radial dimension, r PSD . We set the values of the following parameters to: Furthermore, we assume that the 3000 glutamate molecules are released as a disk of uniform concentration with radius r=0.026 μm.
Several physical parameters related to receptor function can be computed, including the current density J i (r,t) generated by each open state: Here we assume that V m is a constant. We use the following values for the conductance [17]: The total current density J(t) generated by the receptors is simply the sum of the individual J i (t). Since it is difficult to experimentally resolve the radial dependency of the current density, we also will be Parameter Value considering the current I i (t), i.e., the current density J i (r,t) integrated over the active region for each open state, to be given by: The total current I(t) generated by the receptor is simply the sum of the individual I i (t). We give examples of both J(r,t) and I(t) in Figure  3. The figure shows that the current density is peaked at some finite time at the origin (r=0) and rapidly decreases for large radii and longer times. The total current has similar behavior in that it is maximal at some finite time, which we label as t peak .
Since the current I(t) is generically maximal for the single vesicle release event, a useful quantity to define is the maximum current I peak attained at a time t peak over the time series: The total charge generated by the receptor distribution is given by: with the total charge Q being the sum of the individual charges.
While there is no constraint for the receptor density in the equations we use above, the dimensions of the receptors and of the postsynaptic density impose some boundary conditions for this parameter. Thus, from electron microscopic studies, AMPA receptors are generally assumed to be cylindrical proteins with a diameter of approximately 15 nm [26], and the maximum receptor density that can be physically achieved in a postsynaptic density is about 5.7×10 −3 nm −2 : Therefore, this value is highlighted in all our results and graphs (in the form of a dashed line). Furthermore, although this physical value is fixed, it is important to keep in mind that the dynamics of the system depend only on the ratio of N/h D . Therefore, if we double our choice of the Debye length, we would see the same dynamics as our original Debye length with twice the receptor number. Note that in this study, we ignored the phenomenon of glutamate uptake (by astrocytic or neuronal glutamate transporters [27][28][29]), which significantly affects glutamate concentration in the synaptic cleft. However, incorporating this phenomenon would have a similar impact for both the free diffusion and the receptor-limited diffusion case (in fact, it would further emphasize the differences we observe, as shown below).

Free diffusion
In order to understand the differences introduced by the receptor binding/unbinding, we first consider the case in which the receptors do not modify the free glutamate concentration. This amounts to setting F (y) to zero in equation (9). We refer to this as the free diffusion case since glutamate diffuses as if no receptors are present. The free diffusion case can be understood as the limit of N→0 of the rate equations (10). There cannot be any current generated in the free diffusion case since the number of receptors is technically zero. However, one way to proceed is to simply consider equations (10) with N=0 but then use the definitions of J i (r,t) and I i (t) with a finite value of N. This can be understood as decoupling the receptors from glutamate dynamics.
Under these conditions, the behavior of I(t) and Q are shown in Figure 4. The linear relationship between current and charge and receptor density is simply due to the fact that the current is proportional to the number of receptors activated. For a fixed density of receptors, as radial dimension increases, the number of receptors that contribute to the current increases, thereby increasing the rate of change of current (and charge), which accounts for the increasing slope we observe. However, for sufficiently large radial dimension, the receptors located far from the release site experience a small enough glutamate concentration that the open states are barely populated. Therefore, the slope value should eventually asymptote to some maximum slope value.
It is also interesting to consider the individual contribution of the various open states to the total current and charge ( Figure 5). The key point is that, depending on the radial dimension, the open state that contributes most to the total peak current varies. For example, for a 42 nm radial dimension (a relatively small radial dimension), the O3 state contributes most, and then the O4 state, and finally the O2 state, independently of receptor density over the range we study. However, at 162 nm, the O2 state contributes the most, then the O3 state, and finally the O4 state. This is of course not surprising since for sufficiently large radial dimension, most of the receptors lie further away from the center, where glutamate concentration is lower, and the majority of the receptors are using the O2 state to generate current, a state in which binding of only two molecules of glutamate is sufficient to induce channel opening. A similar behavior occurs for the individual open state contribution to the charge, where for small radial dimension, the O3 state makes the largest contribution, while at much larger radial dimensions, and the O2 state has the largest contribution, independently of receptor density. Furthermore, we observe that the behavior of I O2 and Q O3 for the different radial dimensions appears to vary significantly. This apparently strange behavior will be clarified when we include the effects of binding/unbinding on glutamate concentration and receptor dynamics.

Receptor-limited diffusion
We then incorporate the receptor-limited diffusion to perform the same calculations as for the free diffusion case. An example of a direct comparison between the free diffusion and receptor-limited diffusion for peak current and total charge is shown in Figure 6. Interestingly, even for a relatively small number of receptors (<50), there is already a significant difference between the two results. In this example we use a postsynaptic density with a radius of 102 nm, and for N=20 there is already about a 30% difference (corresponding to a difference of 1.3 pA for V m =−70 mV) between the currents generated with the free diffusion and the receptor-limited diffusion case. Figure 7 shows the peak total current and the total charge for different radial dimensions. For very small receptor density, we observe the same linear relationship as in the free diffusion case, which is expected since the limit of small density corresponds to the free diffusion case. However, with increasing receptor densities, we observe very strong deviations from the behavior obtained with the free diffusion case. For all radial dimensions studied, both peak current and charge reach a maximum; although at somewhat different values of receptor density (the charge reaches a maximum at a higher receptor density). As receptor density increases, I peak and Q reach a maximum value and then de-crease. For smaller radial dimensions, this maximum is reached at many higher receptors densities than for larger radial dimensions. This is not surprising since for small radial dimensions, receptors are exposed to a very high glutamate concentration, and thus a very high receptor density is required to observe a receptor-limited diffusion effect. This also explains why the largest I peak and Q values are observed for the smallest radial dimension. This maximum value for I peak and Q at a given receptor density represents therefore the optimal receptor density value for a given radial dimension.
Considering the critical physical receptor density, it appears that the only receptor distributions that can generate a maximum peak current are those with radial dimensions approximately greater than 90 nm, while those that can generate a maximum charge are those with radial dimensions approximately greater than 100 nm. These values correspond to 145 receptors for the former (90 nm) and 180 receptors for the latter (100 nm). (Recall the discussion at the end of section 4 about the close relationship of this value with the Debye length.) Furthermore, for large radial dimensions, the optimal receptor density decreases. This means that the optimal configuration of AMPA receptors for large radial dimensions will deviate significantly from a tight packing configuration. This is not unexpected since for larger radial configurations, glutamate concentration far away from the release site is small, and a decrease in receptor density will reduce competition between receptors for the limited number of glutamate molecules.  Again it is interesting to consider the individual contribution of the various open states to peak current and total charge (Figure 8). The results are strikingly different from those shown in Figure 5. One very important difference is that the O2 state can contribute more than the O3 state. This is an important result, since the O2 state is the only state that shows an increasing contribution with increasing receptor density (up to a point). This is perhaps not entirely surprising since any increase in receptor density would mean an immediate competition between receptors, which would result in a decline in O3 and O4 states' contributions, and allowing the O2 state to play a bigger role. It is exactly this interplay that results in extreme values for peak current and charge.
It is also interesting to analyze the unique patterns observed for I O2 and Q O2 . For large radial dimensions, the maximum of the curve occurs very close to the zero density of receptors (for very large radial dimensions, there is not even a true maximum), and this result accounts for the curious behavior we previously discussed for the free diffusion case. For example, for r PSD =142 nm, this maximum occurs very close to zero, while for r PSD =162 nm, there is no true maximum, and therefore Receptor-limited diffusion case: Effects of receptor density on peak current and total charge. Peak current (A) and total charge (B) were calculated as a function of receptor density. The graphs illustrate the relationship for the following radial dimensions: r PSD =42 nm (purple), r PSD =62 nm (blue), r PSD =82 nm (cyan), r PSD =102 nm (green), r PSD =122 nm (yellow), r PSD =142 nm (orange), and r PSD =162 nm (red). The dashed line corresponds to the maximal receptor density for the AMPA receptor.
at zero density Q O2 /N is greater for the receptor configuration with r PSD =162 nm.

Comparison to experimental results
We compare the results obtained with our two models to published experimental data involving one or two release events ( Figure 1A) [30]. Experimental results show that for a single release event, the peak current (at a potential of −75 mV) is 5 pA, and the peak current doubles to 10 pA for two simultaneous release events. Unfortunately, available experimental data do not provide the number of glutamate molecules released, the density of AMPA receptors, and the radial extent of AMPA receptors. Thus, we choose an AMPA receptor configuration that allows us to fit the experimental data with both the free diffusion and receptor-limited model. We use r PSD =102 nm for both models, and the number of AMPA receptors N=28 for the free diffusion model and N=40 for the receptor-limited model. We present the results for one vesicle release in Figure 9A, where it is clear that the current responses from both theoretical models are almost identical, although they have faster rise and decay times than the superimposed experimental results from reference [30]. All parameters for our AMPA receptor model and diffusion algorithm were optimized [19] with respect to experimental results from reference [23], and the differences observed in Figure 9A are likely due to limitations in experimental recording and dendritic integration, whereby membrane capacitance slows the response. Furthermore, the kinetic rate constants have not been optimized to fit the experimental data, which can dramatically affect the time series. However, the key point here is that all the curves exhibit the same peak current. Using the same AMPA receptor configuration, we calculated the responses with twice the amount of glutamate release to mimic the simultaneous release of two vesicles. Only the receptor-limited model reproduces the doubling of the peak current ( Figure 9B), which is in agreement with our interpretation that receptors in the receptorlimited model are less saturated and therefore able to bind more glutamate molecules.
Although these results suggest that the receptor-limited model provides a better match with experimental data, this finding may not be true for all AMPA receptor configurations. For example, if for a particular configuration, both models result in low saturation of receptors after a single release event, it would be expected that both models would provide similar responses for two release events. Therefore, an important way to differentiate the two models is by controlling the number of receptors, i.e., changing the number or density of receptors while holding the number of glutamate molecules released fixed. As shown in Figure 6, our simulation results predict a substantial difference in the observed peak current as the density of receptors is scaled.

Discussion
Our results indicate that taking into account the feedback controlled receptor-limited diffusion to calculate free glutamate concentration in the cleft significantly changes the patterns of important functional parameters of receptor function, such as total current and total charge through the receptor channel. The physical dimensions at which these effects becomes relevant are within the values experimentally observed (for a reasonable choice of Debye length), indicating that such consideration needs to be incorporated when modeling neurotransmitter/receptor interactions. Furthermore, our results indicate that considering the effect of binding/unbinding shed new lights on the understanding of receptor density and receptor distribution profile within the postsynaptic density. In particular, we discover that, for a uniform distribution of AMPA receptors on a disk, and for a given radial dimension, there is an optimal radial density profile of receptors where peak current or total charge generated by the receptors is maximal, which is a property that cannot be observed using a free diffusion model. This is a direct result of the competition between receptors for the limited number of glutamate molecules, and as expected, when the dimension of the postsynaptic density increases, the optimal receptor density decreases.
Using a typical size for the AMPA receptor, we determine a maximum receptor density limit, and our results indicate that for a given receptor density below this limit, there is a corresponding radial dimension for which peak current or peak charge is maximal. These values are about 90 nm for the former and 100 nm for the latter. While these values are reasonably close to the value of a typical postsynaptic density (between 50 and 500 nm), they are clearly on the low side compared to experiments. However, it is important to note that we have only considered the presence of AMPA receptors in the postsynaptic density and inclusion of NMDA receptors (which have higher affinity for glutamate than AMPA receptors) as well as glutamate transporters will modify this optimum value. The density of receptors at which the maximum occurs corresponds to approximately 145 or 180 receptors (for a Debye length of 1 nm) depending on whether we consider the maximum of peak current or charge.
We made a number of simplifying assumptions throughout this study. First, we assumed that the dynamics of glutamate diffusion perpendicular to the plane (i.e., across the cleft) plays a negligible role. This means that we assumed an instantaneous travel for glutamate upon release across the cleft. In previous studies, it was shown that this effect is negligible [22], and this assumption also greatly simplifies the numerical calculation. However, including these effects would make the role of the Debye length far more realistic, since the receptor would now experience a concentration gradient transverse to the plane. Second, the smearing of the receptor was a necessary assumption in order to solve the problem in the continuum limit. We tried to take into account the physical size of the receptor by focusing our analysis on realistic receptor densities, but it would be interesting to construct models that include finite-size corrections, while retaining the simplicity of the continuum limit. Third, the glutamate-receptor interaction is strictly local, meaning that a receptor at position x (in the plane) can only interact with glutamate molecules (therefore concentration) at position x. However, the finite Debye length implies that the receptor can interact with glutamate concentration at x ± h D as well. It would be interesting to study the effect of such non-local interactions in these models.
An important step forward will be to study a more realistic configuration for the receptors. It is generally assumed that AMPA receptors are distributed in a mosaic pattern at the majority of excitatory synapses [31]. The approach we have used in this study is general enough to accommodate such a distribution, and therefore studying the system integrating NMDA receptors, AMPA receptors, as well as metabotropic receptors and glutamate transporters is an important future direction for our work. Since AMPA receptors in this more realistic situation will be exposed to even less glutamate molecules, the effect of receptor binding/unbinding should be quite relevant. Furthermore, our results indicate that the optimum receptor distribution is likely not uniform. From our results, it is more likely that a tightly packed configuration with AMPA receptors close to the release site and decreasing density away from the release site might correspond to the optimum distribution. This result appears to agree with experimental observations for receptor distributions [32]. It would be very interesting to demonstrate that the integration of the receptorlimited diffusion provides a physical solution for the actual optimal receptor distribution in the postsynaptic density and to evaluate how general this finding could be in different neurotransmitter/receptor systems at different synapses. Furthermore, it would be interesting to adapt our techniques to better understand the clustering of receptors observed experimentally [32]. Finally, these results could also account for the generally admitted view that synaptic plasticity is the result of movements of AMPA receptors closer or further from the center of the postsynaptic density [33].