Synaptic bouton properties are tuned to best fit the prevailing firing pattern

The morphology of presynaptic specializations can vary greatly ranging from classical single-release-site boutons in the central nervous system to boutons of various sizes harboring multiple vesicle release sites. Multi-release-site boutons can be found in several neural contexts, for example at the neuromuscular junction (NMJ) of body wall muscles of Drosophila larvae. These NMJs are built by two motor neurons forming two types of glutamatergic multi-release-site boutons with two typical diameters. However, it is unknown why these distinct nerve terminal configurations are used on the same postsynaptic muscle fiber. To systematically dissect the biophysical properties of these boutons we developed a full three-dimensional model of such boutons, their release sites and transmitter-harboring vesicles and analyzed the local vesicle dynamics of various configurations during stimulation. Here we show that the rate of transmission of a bouton is primarily limited by diffusion-based vesicle movements and that the probability of vesicle release and the size of a bouton affect bouton-performance in distinct temporal domains allowing for an optimal transmission of the neural signals at different time scales. A comparison of our in silico simulations with in vivo recordings of the natural motor pattern of both neurons revealed that the bouton properties resemble a well-tuned cooperation of the parameters release probability and bouton size, enabling a reliable transmission of the prevailing firing-pattern at diffusion-limited boutons. Our findings indicate that the prevailing firing-pattern of a neuron may determine the physiological and morphological parameters required for its synaptic terminals.

The morphology of presynaptic specializations can vary greatly ranging from classical single-release-site boutons in the central nervous system to boutons of various sizes harboring multiple vesicle release sites. Multi-release-site boutons can be found in several neural contexts, for example at the neuromuscular junction (NMJ) of body wall muscles of Drosophila larvae. These NMJs are built by two motor neurons forming two types of glutamatergic multi-release-site boutons with two typical diameters. However, it is unknown why these distinct nerve terminal configurations are used on the same postsynaptic muscle fiber. To systematically dissect the biophysical properties of these boutons we developed a full three-dimensional model of such boutons, their release sites and transmitter-harboring vesicles and analyzed the local vesicle dynamics of various configurations during stimulation. Here we show that the rate of transmission of a bouton is primarily limited by diffusion-based vesicle movements and that the probability of vesicle release and the size of a bouton affect bouton-performance in distinct temporal domains allowing for an optimal transmission of the neural signals at different time scales. A comparison of our in silico simulations with in vivo recordings of the natural motor pattern of both neurons revealed that the bouton properties resemble a well-tuned cooperation of the parameters release probability and bouton size, enabling a reliable transmission of the prevailing firing-pattern at diffusion-limited boutons. Our findings indicate that the prevailing firing-pattern of a neuron may determine the physiological and morphological parameters required for its synaptic terminals.

INTRODUCTION
Chemical synapses are specialized compartments of neurons for neuron-to-neuron or neuron-to-muscle communication that exist in various morphological layouts (Rollenhagen and Lubke, 2006) ranging from classical single-release site synapses at cortical neurons (Gulyas, 1993) over intermediate-sized multi-releasesite boutons of the calyx of Heldt (Satzler et al., 2002) or the arthropod neuromuscular junctions (NMJs) (Bradacs et al., 1997) to giant multi-release-site terminals at vertebrate NMJs (Dreyer et al., 1973). While the principal processes underlying synaptic transmission and plasticity at these synapses are rather well understood, the role that the various morphological layouts may play in shaping synaptic functions is currently unclear. One reason for this astonishing discrepancy may be the difficulties to systematically control the morphological and physiological configuration of a given synaptic system.
Prominent examples for two highly related but morphologically distinct multi-release-site boutons are the well-characterized glutamatergic synapses at larval NMJs of Drosophila (Jan and Jan, 1976;Atwood et al., 1993;Schuster, 2006;Collins and DiAntonio, 2007). Larval body wall muscles of Drosophila are typically innervated by two motor neurons of which one forms large spherical and somewhat variable type Ib boutons and a second motor neuron forms smaller and more regular type Is boutons. Both types of boutons harbor multiple glutamatergic release sites that are spaced according to a strict nearest-neighborhood relationship (Atwood et al., 1993;Meinertzhagen et al., 1998;Sigrist et al., 2003) and that show discrete differences in their release probabilities (Kurdyak et al., 1994;Cooper et al., 1995). In spite of their close functional relationship, it has so far not been possible to systematically assess why both motor neurons maintain their obvious morphological and physiological differences and what their relevance is for the animal's behavior.
In order to shed light onto the poorly understood structurefunction relationship of such discrete bouton-configurations, we used a combined theoretical and experimental approach to dissect the parameters that are biophysically important to support the natural function of discrete bouton configurations. We therefore developed and verified a three-dimensional model of synaptic bouton functions of larval NMJs of Drosophila and used it to systematically assess the meaning of experimentally rather inaccessible bouton parameters like bouton size, vesicle diffusion and vesicle release probability. We found that these parameters affect bouton output and its reliability of transmission on different time scales and that their individual combination in real boutons is predictive of the prevailing firing pattern that a given bouton type experiences.

MORPHOLOGICAL AND PHYSIOLOGICAL PARAMETERS TO GENERATE A MODEL OF GLUTAMATERGIC BOUTONS
To systematically assess the influence of morphological and physiological constraints on the performance of multi-release-site boutons we first established a functional reconstruction of such boutons in silico that aimed at simulating the vesicle dynamics within a bouton during trains of action potentials. From confocal image stacks of NMJs of third instar Drosophila larvae ( Figure 1A) and representative measurements of bouton dimensions of type Is and type Ib synaptic terminals, we first estimated the typical bouton diameter of type Is and type Ib boutons to 2 and 3 µm, respectively (right panel in Figure 1A). In addition, some terminal boutons showed diameters of up to 5 µm. These different bouton morphologies were constructed in ProMesh (Reiter and Wittum, 2014), and discretized with a tetrahedral volume grid. The geometries represented spheres with a volume spared out at the center of each sphere representing organelle occupation, hence an area not occupied by vesicles (Figures 1C,D, 2). Each bouton geometry consisted of a zone 0 where vesicles were motile as well as exocytosis zones 1 . . . n where vesicles are motile and could be released ( Figures 1B-D, 2). Taking into account the total number of vesicles that can be released from a NMJ during vesicle depletion experiments (Delgado et al., 2000) and data from ultrastructural analyses (Sigrist et al., 2002) we estimated a concentration of 150-500 vesicles per µm 3 within a bouton at rest (see Materials and Methods for details). The positioning of the vesicle release sites on the surface of these spherical boutons was guided by serial ultrastructural reconstructions of larval boutons from which we determined the densities of T-bar harboring active zones between 0.37 and 0.46 µm 2 (Meinertzhagen et al., 1998;Sigrist et al., 2002Sigrist et al., , 2003 with an average active zone diameter of 300-400 nm (Meinertzhagen et al., 1998;Sigrist et al., 2002Sigrist et al., , 2003. Based on these data we assumed 7 equally spaced active zones per Is bouton, 10 release sites per Ib bouton and up to 14 release sites in a large terminal bouton. Synaptic active zones were represented by a multiedgesurface approximating a cylindrical structure with a diameter of 300-400 nm and a volumetric depth of 200 nm, which defines a volume that harbored roughly 6-8 vesicles, depending on the starting vesicle densities, which were chosen between 150 and  (muscle 13) showing the synaptic terminals of two motor neurons that form type Is or Ib boutons, respectively. Release sites (green) were labeled with an antibody recognizing the active zone protein Bruchpilot (Wagh et al., 2006). The axonal surface (red) was labeled by anti-HRP immunoreactivity. Right panel: magnification of representative Is and Ib boutons (arrow and arrowhead) displaying an average diameter of 2 and 3 µm, respectively. Scale bars: 5 µm (left) 1.5 µm (right). (B-D) Release sites (green) and the vesicular diffusion area (red) are represented in silico by tetrahedral volume grids (Bastian et al., 2000). Release sites are positioned according to nearest neighborhood laws derived from ultrastructural data (Meinertzhagen et al., 1998;Sigrist et al., 2003). The bouton center is free of vesicles (Denker et al., 2009). (E,F) Simulation of vesicle dynamics in a Ib bouton (3 µm) during 30 Hz or 60 Hz stimulation (P o = 7%) with initial vesicle density of 275 µm −3 . Encoded in color is the spatially resolved vesicle density (in vesicles/µm 3 ) at the indicated stimulus number (599-601). Note that 60 Hz-stimulation rapidly depletes most release sites and their immediate surroundings of releasable vesicles (white) while more distant areas maintain high vesicle densities (red). 500 vesicles per µm 3 . An overview of the parameter ranges used in the simulations is given in Table 1. Where simulation data and results are shown, we specify which values within the given ranges were chosen. Within the active zones we made a transition from vesicle motion at the continuous scale to discrete vesicles that undergo exocytosis (see Materials and Methods for a detailed description). At this stage of the model we did not incorporate voltage-gated calcium channels on the plasma membrane and presynaptic calcium dynamics and its role in triggering vesicle fusion but instead modeled action potential triggered release of vesicles by the term f syn that described a stochastic vesicle release process controlled by the stimulation frequency ω and To implement into our model a meaningful simulation of vesicular motilities during trains of high-frequency action potential firing we made the following considerations: at Drosophila NMJs presynaptic vesicles have been assigned to three functionally distinct vesicle pools, the readily releasable pool, the recycling pool and the reserve pool of vesicles (Kuromi and Kidokoro, 2000;Denker et al., 2009). Although the number of vesicles belonging to each vesicle pool and their respective motilities are very different at rest (Kuromi and Kidokoro, 2000), stimulation frequencies above 10 Hz have been shown to be sufficient to recruit vesicles from all pools for release (Kuromi and Kidokoro, 2000;Wagh et al., 2006;Denker et al., 2009). Since the natural firing pattern of larval motor neurons uses typically frequencies higher than 10 Hz (Figure 6) all vesicles within a bouton have to be considered mobile during natural activity patterns in vivo. It is, however, currently not clear whether vesicle motility during high frequency stimulation follows a simple diffusion process (Gaffield et al., 2006) or whether it underlies active transport mechanisms requiring cytoskeletal elements (Ryan, 1999;Jordan et al., 2005;Tokuoka and Goda, 2006;Seabrooke et al., 2010;Kisiel et al., 2011). We have previously shown that the presynaptic actin-based cytoskeleton within Drosophila boutons is highly dynamic with major rearrangements occurring on a rapid time scale (Steinert et al., 2006). These cytoskeletal rearrangements may result in a constant intermixing of presynaptic vesicles and hence in a dynamic and overall rather homogeneous distribution of vesicles within the presynaptic terminal. Based on these and the above observations we first assumed that all presynaptic vesicles are mobile during stimulation frequencies higher than 10 Hz (Kuromi and Kidokoro, 2000). We, secondly, assumed that during such stimulation periods vesicle motility may be approximated by undirected diffusion-like movements (Shtrahman et al., 2005;Darcy et al., 2006;Fernandez-Alfonso and Ryan, 2008;Westphal et al., 2008). We therefore used an effective diffusion coefficient of D = 5 · 10 −3 µm 2 /s (Shtrahman et al., 2005) that included stick-and-diffuse processes that have been shown to occur during interactions of vesicles with the cytoskelleton (Darcy et al., 2006;Fernandez-Alfonso and Ryan, 2008;Westphal et al., 2008).
The above considerations formed the basis for the formulation of the diffusion equation and the sink term that were solved on complex morphologies in three-dimensional space using numerical discretization methods and solvers (see Materials and Methods section for details). All numerical methods were part of the simulation platform uG (Bastian et al., 2000). Numerical simulations of our model generated three-dimensional density profiles of presynaptic vesicles that allowed us to monitor the local distribution and dynamics of mature vesicles as a function of bouton diameter, the applied stimulation pattern and the vesicle release probability (Figures 1E,F, Supplemental Video S1).

DISPROPORTIONAL DEPLETION OF VESICLE RELEASE AT HIGH STIMULATION FREQUENCIES IS DUE TO LIMITED VESICLE DIFFUSION
To verify the validity of the parameters of our bouton model we wished to compare its behavior in a task that is similar to a simple in vivo experiment. We therefore chose an experimental vesicle depletion paradigm that is triggered by continuous highfrequency stimulation of larval motor nerves in the presence of a blocker of the vesicular proton pump, Bafilomycin A1 (Kuromi and Kidokoro, 2000) (BafA1 in Figure 3A). Synapses in BafA1treated larval filet preparations fail to fill recycling presynaptic vesicles with the neurotransmitter glutamate and hence can be used in intracellular recording experiments of postsynaptic membrane potentials to exclusively visualize the release of preexisting mature vesicles. As expected, the amplitudes of evoked excitatory junctional potentials (eEJPs) of BafA1-treated larval preparations depleted significantly faster during continuous 30 Hz nerve stimulation than in untreated preparations (Figures 3B-D). The here observed rate of eEJP-deletion in BafA1-treated preparations was similar to the rate of vesicle depletion in temperature-sensitive mutants of dynamin (shibire ts ) (Siddiqi and Benzer, 1976) that at restrictive temperatures inhibit all endocytic processes (Delgado et al., 2000) indicating that in spite of the different molecular mechanisms of both treatments the measurable outcome was comparable. In addition, we observed that the eEJP-depletion progressed at rather slow rates during stimulation frequencies between 15 and 30 Hz, whereas at higher stimulation frequencies (60-80 Hz) the initial depletion was disproportionately faster (arrows in Figure 3E). This unexpected phenomenon was unlikely due to desensitization of postsynaptic glutamate receptors, since the quantal amplitudes (miniature excitatory junctional potentials, mEJPs) were unaltered under all tested conditions (data not shown). Instead, the latter observation suggested that a presynaptic limitation was responsible for the disproportionally rapid decay of eEJP amplitudes at high stimulation frequencies.
An analysis of NMJs that transmitted naturally produced motor patterns in the absence of BafA1 (Supplemental Figures S1A,D) revealed that these synapses were capable of maintaining their eEJP amplitudes throughout entire patterns (Supplemental Figures S1C,F) even when the firing frequencies were in the range of 60-80 Hz (Supplemental Figures S1B,E). This result demonstrated that high-frequency firing patterns could be faithfully transmitted by the glutamatergic boutons of larval NMJs. It further suggested that the disproportionally rapid decay of eEJP amplitudes seen in BafA1-treated synapses might be due to limited availability of preexisting mature vesicles at active zones. If such a limited availability of transmitter filled vesicles would cause the disproportional synaptic depression during 60-80 Hz stimulation rates one would expect that short pausing intervals between stimulation trains allow distant mature vesicles to functionally reoccupy release sites over time resulting in recovered eEJP amplitudes. If BafA1-treatment would potentially unspecifically block other synaptically relevant processes such as calcium export or cytoskeletal organization that would disproportionally stronger affect high frequency transmission, eEJP amplitudes should continue to show depression after a paused stimulation. This assumption is contradicted by our results in Figures 3F,G showing transient depression at all stimulation frequencies excluding a persistent pleiotropic defect induced by the procedure. In experiments in which trains of 100 stimuli were separated by pausing intervals of 5 s we found an almost complete recovery of the eEJP amplitudes at all used FIGURE 3 | Disproportional depletion of eEJPs is due to limited availability of presynaptic vesicles. (A) Experimental design: larval filet preparations were pre-incubated in normal saline (NS) or NS containing 2 µM Bafilomycin A1 (BafA1), a blocker of the vesicular proton pump, to inhibit transmitter filling in recycling vesicles. Motor nerves were then continuously stimulated with a suction electrode at various frequencies (HFS) while we simultaneously recorded the membrane potential of the postsynaptic muscle fiber 6. (B,C) Representative traces of these recordings in the absence or presence of BafA1 taken at the indicated time-points of stimulation (stars: stimulation artifacts). (D) Quantification of eEJP-depletion dynamics during continuous 30 Hz stimulation with and without BafA1 incubation. Note, that eEJPs could be detected over at least 15 min of 30 Hz-stimulation in controls without BafA1-treatment (B and circles) whereas eEJPs of BafA1-treated preparations depleted within 2 min (C and squares). (E-G) eEJP-depletion dynamics in BafA1-treated preparations caused by the indicated stimulation frequencies and expressed as a function of the stimulus number. (E) Shows the eEJP-depletion dynamics during continuous long-term stimulation. In (F,G) stimulation trains of 100 stimuli were separated by pausing intervals of 5 s. Note that the depletion of eEJP-amplitudes was disproportionally faster during 60 and 80 Hz stimulations compared to 15 and 30 Hz stimulations. Short pausing intervals resulted in an almost complete recovery of eEJP-amplitudes at the beginning of the experiment (F). In later phases of the experiment (G) the eEJP recovery was significantly higher in 80 Hz stimulated NMJs than in 15 Hz stimulated NMJs (p<0.033 at stimulus 1800 and p<0.019 at stimulus 1900). Data in (D-G) represent means ± SEM of 4 to 6 independent preparations. stimulation frequencies (15-80 Hz) immediately after the pause ( Figure 3F). Furthermore, the longer this mode of stimulation continued the more the recovered eEJP amplitudes deviated between high and low frequencies ( Figure 3G): eEJP recovery was significantly larger in boutons after 1800 stimuli at 80 Hz than at 15 Hz (p < 0.033, n = 4) or after 1900 stimuli (p < 0.019, n = 4). Both observations were consistent with the idea that vesicle motility was limiting the availability of mature vesicles at active zones. This leads during high frequency stimulation to a rapid failure to reoccupy empty active zones with mature vesicles and thus to a rather small total number of released vesicles per train. The higher number of remaining mature vesicles is then responsible for the strong eEJP recovery even after prolonged stimulation periods (triangle in Figure 3G). Conversely, during lower stimulation frequencies vesicle release and active zone reoccupation appear to be more balanced resulting in a higher total vesicle output per train and consequently in a weaker eEJP recovery after prolonged stimulation (circle in Figure 3G).

Verification of the model
To further verify the quality of our model we systematically compared simulated and experimental vesicle depletion curves of different conditions (Figure 4). For the model we used as before parameter values from previous experimental studies (Kurdyak et al., 1994) and set up our model with 6 Ib and 13 Is boutons with P o values of 7 and 29% for the release sites of type Ib and type Is boutons, respectively (Kurdyak et al., 1994). Under these conditions the in silico simulations of the vesicle depletion experiments approximated the experimental observations in that the

FIGURE 4 | Model validation at different frequencies and stimulation protocols. (A)
Comparison between experimental and simulated data at the frequencies 15, 30, 60, and 80 Hz. In order to compare simulation and experiment, the initial vesicle density was computed from the experimentally recorded vesicle release for each frequency and thus initial vesicle densities were set to 500, 280, 190, and 150 vesicles/µm 3 respectively in the simulations. Note that in the beginning 1000 stimuli at 15 Hz the depletion is considerably faster in real boutons than in the simulation. This may indicate that 15 Hz stimulation is initially not sufficient to efficiently mobilize all presynaptic vesicles e.g., from the reserve pool of vesicles. This limitation is overcome at higher stimulation frequencies. (B-E) Experimental (dotted) and simulation (lines) data in stop and go stimulation protocol, with 100 stimuli and 5 s pause, show decline in in the postsynaptic amplitudes (B,C: simulation data averaged over 25 stimuli; D,E: simulation data averaged over 5 stimuli). The plots show that with increasing frequency (from B-E) amplitudes decrease slower, due to higher exocytosis failure rates at high frequencies caused by diffusion limited active zone replenishment. Note that similar to (A) the experimental depletion is faster than the in silico depletion. As above this may be due to a not efficient mobilization of presynaptic vesicles at 15 Hz.

Frontiers in Computational Neuroscience
www.frontiersin.org September 2014 | Volume 8 | Article 101 | 5 simulated total amount of vesicular release decreased at various stimulation frequencies with similar time courses as in the in vivo experiments (Figure 4). Intriguingly, we identified the same two features that together could explain the disproportionally faster depletion of vesicle release at high stimulation frequencies with our model simulation of the local vesicle dynamics within a bouton ( Figures 1E,F): first, 600 stimuli applied at 60 Hz depleted most release sites and their immediate surroundings of releasable vesicles (white patches in Figure 1F) whereas the same number of 30 Hz-stimuli left more release sites functional with one or more releasable vesicles ( Figure 1E). Second, 600 stimuli at 60 Hz caused a decrease of the overall number of vesicles per bouton to about 84% (mostly red areas in bouton lumen of Figure 1F), whereas the same number of stimuli applied at 30 Hz resulted in a stronger overall decrease to about 73% (mostly yellow in Figure 1E). Since we modeled vesicle dynamics during synaptic stimulation as an undirected diffusion process that matched our experimental observations remarkably well, it seemed likely that vesicle diffusion is an important limiting factor during high frequency stimulation episodes in vivo (see Discussion). This included the disproportionally faster depletion of vesicle release at 60 or 80 Hz compared to lower frequencies ( Figure 4A). We further compared experimental and simulation data in stop and go experiments (from Figures 3F,G), where we applied 100 stimuli at the given frequency and 5 s stimulation pause alternatingly (see also Figures 3F,G). Data shown in Figures 4B-E showed that the peak amplitudes decreased with increasing stimulus number, as expected for all stimulation frequencies. We found that both experiment and simulation showed a rapid recovery from depletion during the pausing periods (Figures 4B-E). In both cases the recovery was considerably larger during high-frequency stimulation (Figures 4D,E) than at lower frequencies ( Figures 4B,C). This could be attributed to the fact that the overall failure rate in vesicle release (i.e., a release event is triggered, but no vesicle is present at the active zone) was larger at high frequencies, which in response causes a slower decline in postsynaptic response. This effect was a further indicator for vesicle diffusion limiting synaptic response. In order to quantify the results in Figure 4A, we applied different measures (least-squares and L 2 -norm, see Materials and Methods for more detail) to the experimental and simulated data ( Table 2). These similarities demonstrate that our model includes the necessary biophysical laws that explained the stimulationinduced vesicular dynamics of multi-release-site boutons in space and time and hence might be used to analyze vesicle movement and their potential restriction for release at high temporal and spatial resolution. We saw that the overall bouton dynamics in vivo and in silico show significant similarities, though the diffusion-based in silico results demonstrated a slower decay phase at high frequency stimulation in the initial stimulation phase, thus indicating that pure diffusion might not suffice for high amplitude/high frequency output (see Discussion).

Bouton performance is differentially affected by bouton size, release probability, and stimulation frequency
Since our three-dimensional model described the vesicle dynamics in a multi-release-site bouton during example stimulation  We found that high P o -values strongly affected the magnitude of the overall vesicle release at the beginning of a 40 Hz stimulation train ( Figure 5A). However, the strong initial vesicle release disappeared within a second of stimulation due to the rapid depletion of vesicles at individual release sites. Thus, high P o 's seemed only meaningful if large amplitude outputs were needed for short time periods, whereas small P o 's provided long-term reliability of smaller amplitude outputs ( Figure 5A). The size of boutons affected its performance under otherwise constant conditions (P o = 6%; ω = 40 Hz) on a rather long time-scale. Bigger boutons with correspondingly more release sites showed larger outputs that depleted significantly slower than those of smaller boutons ( Figure 5B). Thus, larger boutons seemed to be better suited for sustained periods of synaptic activity whereas smaller boutons with fewer release sites were only reliable during short trains. Note that the term "reliable transmission" refers to the phases in which vesicle exocytosis is constant, thus the phase in which the system is not limited by high stimulation frequencies, bouton size or vesicle output probability.
Our analysis of the effects of a broad range of stimulation frequencies revealed that frequencies of up to 40-50 Hz were reliably transmitted by our standard type Ib bouton for seconds (initial plateaus in Figure 5C). At higher stimulation frequencies reliable vesicle release was limited to only few stimuli until it becomes diffusion-limited (see Figure 5C). Thus, type Ib boutons could transmit a broad range of stimulus frequencies even in the absence of vesicle recycling, however, frequencies above ≈50 Hz could only be maintained for very short time-periods ( Figure 5C).
The data in Figure 5B, with initial vesicle density set to 400 µm −3 and 7 active zones for the Is bouton, show that the exocytosis dynamics are separated into a stable phase in which exocytosis is constant and an unstable phase in which there is a decay in vesicle exocytosis. The law in Figure 5C shows that in the range above 50 Hz the system becomes unstable nearly instantaneously. Experimentally this effect is observed at around 45 Hz for an output probability estimated at roughly 7%. FIGURE 5 | Bouton output and output endurance are differentially affected by stimulation frequency, bouton size, and P o . Vesicle density is set to 400 µm −3 . In each figure row only one of these three parameters were varied, while keeping the other two fixed. Systematic analysis of the effects of changes in P o (A), bouton diameter (d = 2, 3 and 5 µm) with 7, 10, and 14 release sites respectively (B) and stimulation frequency (C) on the simulated total output dynamics over several simulations (dotted) or on derived biophysical depletion laws (lines). If not indicated otherwise we used a P o of 6%, a bouton diameter of 3 µm and a stimulation frequency of 40 Hz. Note that the total bouton output is initially very large at high P o -values, but it rapidly depletes to the output of low-P o boutons that continue reliable transmission. Further note that larger boutons harbor more release sites resulting in larger and longer-lasting total outputs (B). Stimulation frequencies up to 40 Hz can be transmitted reliably over some time (plateaus in C) while higher frequencies can only be maintained for short time periods. Interestingly, the long term behavior is independent on the choice of P o for P o ≥ 5%.

Decomposition of the compound response of Is/Ib-harboring NMJs
Since all three analyzed parameters affected bouton performance in distinct manners it was tempting to combine these individual model-derived characteristics into realistic functional profiles of the two larval bouton types. If one assumed that the active zones of the smaller type Is boutons operated with a rather high vesicle release probability (in the range of P o = 40%) while the bigger type Ib boutons used a lower, in the range of P o = 6%, (Kurdyak et al., 1994;Lnenicka and Keshishian, 2000) our above simulations allowed clear predictions as to how the output of a NMJ with roughly similar numbers of both bouton types will look like ( Figure 6A): a continuous stimulation of, e.g., 40 Hz with initial vesicle density set to 400 µm −3 will result in a compound junctional output (dark gray line) that is initially dominated by the strong release from Is boutons (black line) before their contribution decays due to rapid vesicle depletion at their active zones (white patches in Figure 6B). In contrast, the type Ib boutons of this NMJ are expected to show reliable but small amplitude performance (gray line in Figure 6A) with no extensive signs of diffusion-based vesicle depletion (left bouton in Figure 6B). Intriguingly, in our in vivo experiments in which we stimulated the axons of both motor neurons synchronously at various frequencies ( Figure 3E) we observed a very similar behavior of the compound response as indicated by an initial steep decrease in the eEJP amplitudes followed by an abrupt transition (arrows in Figures 3E-G) into a slower decay phase. These results suggested both in silico and in vivo that type Is boutons were not well suited to reliably transmit sustained high-frequency stimuli and they indicate that it was unlikely that in a behaving animal the natural motor pattern fires type-Is and -Ib boutons synchronously.

The natural firing pattern matches the functional properties of type Is and Ib boutons
Our above results predicted that Ib boutons are likely utilized to transmit longer-lasting high-frequency stimuli whereas Is boutons are better suited for few low-frequency events. We tested these predictions by analyzing the natural firing pattern that is generated by the larval central pattern generator to trigger coordinated muscle contraction. In order to differentiate between the patterns transmitted by Is and Ib boutons we established simultaneous recordings from the neighboring muscles 6 and 13 that shared a common Is-innervation from the same motor neuron but received an independent Ib-innervation each from two distinct motor neurons (Lnenicka and Keshishian, 2000). We found that the natural motor pattern consists of two discrete components ( Figure 6C): a tonic eEJP pattern originating from type-Ib boutons (arrowheads) that displayed in each train a bell-shaped frequency ramp of pulses ranging between 10 and 45 Hz on muscle 6 (black line in Figure 6D; 50.8 ± 5.4 pulses per train,  Initial vesicle density was set to 275 µm −3 . This NMJ exocytosis pattern was driven by the brain governed motor pattern as measured experimentally and presented in (D).
mean ± s.e.m. of n = 10). The second component originated from type-Is boutons and consists of only few large pulses at a rather stable frequency of about 12 Hz (arrows in Figure 6C; gray line in Figure 6D; 5 ± 0, 7 pulses per train, n = 10). These results demonstrated that Is or Ib motor neurons fire with patterns that corresponded well to our above predictions of their individual physiological properties (see Discussion and Figure 7). Conversely, when we fed the naturally observed firing pattern of both boutons into our model we found that the in silico boutons reproduced the respective patterns consistently even in the absence of vesicle recycling ( Figure 6E).

ELECTROPHYSIOLOGICAL RECORDINGS
Electrophysiological recordings of membrane potentials of body wall muscle 6 were performed on filet preparations of sizematched mid third instar male larvae at 22C. Dissected larvae were first immersed in normal saline (NS) solution containing 2 mM Ca 2+ (Jan and Jan, 1976) that was then replaced by NS containing 2 µM Bafilomycin A1 (BafA1). After 5 min incubation muscle 6 of segment A 2−3 was impaled with a 1530 M microelectrode filled with 3 M KCl. For stimulation the cut end of the intersegmental nerve of this segment was placed into a suction electrode and suprathreshold stimuli were applied at the indicated frequencies with or without intermitting pausing intervals. Cells with a resting potential more negative than −60 mV and an input resistance R in of 5 M were used for further analysis. Data analysis was performed off-line (Clampfit10, Axon Instruments). Recordings of natural motor patterns were performed similarly (NS containing 2 mM Ca 2+ without BafA1) except that the motor nerves were left attached to the ventral nerve cord and that we recorded simultaneously from muscle 6 and muscle 12 or 13. This configuration allowed us to differentiate the firing patterns of type 1b-and 1 s-innervations (Chouhan et al., 2010). Only those trains were taken for further analysis that correlated with a peristaltic contraction wave of the segmentally repeated body wall muscles. eEJP-amplitudes in Supplemental Figure S1 were measured from the extrapolated resting potential to the peak of the eEJP.

Bouton size
Paraformaldehyde fixed larval filet preparations were immunolabeled with antibodies recognizing a neuronal surface epitope (anti-HRP) and the active zone protein Bruchpilot (Wagh et al., 2006). Confocal images of 12 of such labeled M12/13-NMJs (e.g., Figure 1A) were used to estimate the diameters of synaptic boutons of type Is motor neurons and type Ib motor neurons using the measurement tool of ImageJ. As exemplified in Figure 1A bouton sizes varied more in type Ib branches than in type Is branches but they were on average larger than type Is bouton. Terminal boutons showed particularly large diameters. Because of these considerable variances we chose to confine our current analysis to three representative classes of bouton sizes with the following diameters: Is: 2 µm; Ib: 3 µm; terminal bouton: 5 µm.

Vesicle density
The vesicle densities used in our model were determined as follows: (1) The total number of quanta that can be released from a Drosophila NMJ has been previously determined to be in the range of 85,000 (Delgado et al., 2000). This translates to 850-1700 vesicles per bouton for NMJs with 50-100 boutons and a concentration of 130-270 vesicles/µm 3 for an average bouton (diameter: 2.5 µm, inner vesicle free diameter: 1.5 µm).

Output probability P o
The output probability P o values used in this study correspond to the experimentally determined vesicle release probabilities P of synapses of the two bouton classes Is and Ib that were estimated to average P values of 29 and 7%, respectively (Kurdyak et al., 1994;Cooper et al., 1995). Although a large body of evidence suggests that vesicle release probabilities are non-uniform and subject to change during facilitation and depression phenomena we chose in this stage of our model to simplify and keep the P o values constant during each simulation. To still account for such changes we assayed the effects of different P o values in a bandwidth of 1-90%. The P o in our model describes the overall likelihood that a given vesicle is released from a given active zone and therefore represents a summed value for all processes involved in an individual vesicle release event including vesicle docking and priming, the functional status of voltage gated calcium channels, presynaptic calcium dynamics and calcium triggered vesicle fusion.

DIFFUSION-REACTION MODEL FOR VESICLE DYNAMICS
The mathematical model presented in this manuscript is based on the coarse-scale model of brownian particle motion proposed in Einstein (1905). Using this approach we can deduce that vesicle dynamics can be described by means of a diffusion-reaction equation. To achieve this result, a fixed control region V ⊂ R 3 and an interval of time I = [0, T] are considered. The total number of vesicles contained in V at time t ∈ I is calculated as where c (units [c] = vesicles · µm −3 ) is the concentration of the vesicles. The variation of N V over time is regulated in principle by the flux of vesicles through the boundary ∂V, a source and a sink term, i.e., Assuming regularity, we can extend J to the interior of the control volume V. Using Gauss' Theorem, and a standard localization argument, yields the local form of (2): In (2) and (3), J represents the vesicle flux vector, while R and A represent, respectively, recycling and absorption of vesicles. In general, both R and A depend on concentration and may also exhibit explicit dependence on time and spatial coordinates. Since A should vanish when concentration is zero, it is given through an expression of the type A(c, x, t) ≥ 0, such that the condition A(0, x, t) = 0 applies. The physical units of A are [A] = vesicles/(µm 3 s). For the purposes of our manuscript, recycling is switched off (i.e., R is set equal to zero from the outset), and only the depletion of vesicles is accounted for. Using Fick's Law, and assuming isotropic diffusion, the vesicle flux vector can be expressed by J = −D∇c, with D the diffusion coefficient, and Equation (3) can be rewritten as The computational domain, over which (4) is solved, is denoted by ⊂ R 3 and coincides with the region of space occupied by a bouton. The set is partitioned in subdomains i , i = 0, . . . , S, with S being the number of synapses. 0 is the subdomain containing vesicles which are not part of synaptic regions and accounts for the largest part of the bouton. For i = 1, 2, . . . S, the i are the synaptic regions. The computational domain therefore is written as see Figure 2. For all outer boundary surfaces, we use constant Neumann zero boundary conditions: J · n = −(D∇c) · n = 0, ∀ x ∈ ∂ , and ∀ t ∈ I. (6) Finally, we close the mathematical problem by prescribing the initial condition where t = 0 denotes the origin of time.

THE SINK TERM
The reason we use a sink term rather than a Neumann flux boundary condition to simulate vesicle exocytosis is that such a boundary condition would describe the situation in which the vesicles are permitted to cross the membrane of the bouton. However, this is not the observed mechanism of vesicle depletion in a bouton. Indeed, when the signal of the action potential reaches a synaptic region, the vesicles available for exocytosis fuse with the inner surface of the bouton and release neurotransmitters, which, in turn, are allowed to cross the membrane of the bouton. Therefore, in spite of a zero Neumann boundary condition for the vesicle flux, the total number of vesicles contained in a bouton is not conserved over time. To account for the fusion of the vesicles with the membrane of the bouton, we introduce a sink term f syn ≡ A that, under the conditions discussed below, removes one or no vesicle per action potential from each synaptic region of the bouton. Vesicle exocytosis is a discrete event that reduces the pool of vesicles by an integer number of exocytosed vesicles. The sink term defined in the following models exocytosis based on the general synaptic vesicle release probability P o and a stochastic probability P i that is calculated for each synapse in each simulation time step. Since vesicle motion is modeled on the continuous scale, and vesicle exocytosis relates to single vesicles, the sink term needs to couple the continuous scale with the single particle scale. This leads us to the following definition of the vesicle releasing sink term f syn : with x ∈ , t n = n ω ∈ I, and n = 1, . . . , N A . In the following, we explain the terms used in the above equations: 1. Given the interval I, that represents the lapse of time during which the system is observed, the natural number N A ∈ N is the total number of times (modeled as points of I) in which the action potential stimulates the depletion of vesicles. The generic instant of time t n ∈ I, with n = 1, . . . , N A , corresponds to the nth stimulus, and is defined by t n := n ω , with ω being the stimulation frequency. 2. i refers to the synapse number. 3. The mathematical definitions of the Dirac Delta δ, Heaviside distributions and + , and characteristic function χ i are given in Section 3.9. The Dirac Delta δ (t − t n ) restricts exocytosis to the time points t n of action potentials; the characteristic function χ i restricts the activity of the sink rate solely to the ith synaptic region. In the present setting, the Dirac Delta and the Heaviside distributions are dimensionless, while h(x, t − n ) has units µm −3 . 4. P o defines the release probability of vesicles at single synapses. 5. N i (t − n ) − 1 is introduced to ensure that the Heaviside function + is equal to one, if there is at least one vesicle present at the corresponding release site. Therefore, in order to obtain the result + N i (t − n ) − 1 = 1, the quantity N i t − n − 1 has to be non-negative. 6. A random generator produces values P i between 0 and 1 for each synapse at every action potential. If P i is lower than P o , a vesicle is exocytosed. This is regulated by the expression Frontiers in Computational Neuroscience www.frontiersin.org September 2014 | Volume 8 | Article 101 | 10 (P o − P i (t n )) summed over all synapses i. Thus, in between action potentials the model describes a standard diffusion process. 7. Given a generic function of time, ϕ : I → R, the notation ϕ t − n stands for ϕ t − n = lim t → t − n ϕ(t).
Moreover, the constant a 0 shall be set equal to unity from here on. The sink term (8) is now treated in the following way. We conventionally set t 0 ≡ 0 and t N A + 1 = T, and write the interval In each open interval (t m − 1 , t m ), m = 1, . . . , N A + 1, f syn (x, t) vanishes identically and, thus, the solution to Equation (4) represents pure diffusion. However, if at a given t p = p ω ∈ I, with p ∈ {1, . . . , N A }, the Heaviside functions featured in (9) are simultaneously different from zero for the synaptic region i , the sink term introduces a jump in the value of concentration at t p , for all x ∈ i . Hence, by employing the notation c( The "updated" concentration c x, t + p is treated as initial condition for the diffusion equation in the interval (t p , t p + 1 ).
We remark that, according to (11), only one vesicle is removed from the region i at t p . Indeed, by integrating (11) over i , and recalling that N i t ± p = i c x, t ± p dv, we obtain We denote now by c n , with n = 1, . . . , N A , the restriction of the solution to Equation (4) to the interval (t n − 1 , t n ). By continuating each c n to the time point t n , we set c n (x, t n ) ≡ c(x, t − n ), and define the number of vesicles contained in i at time t n in the following way: To evaluate the Heaviside function + of Equation (9), we use the number of vesicles N i (t n ), determined in Equation (13). More generally, however, for a given x ∈ , the jump in the vesicle concentration at time t n , with n = 1, . . . , N A , is given by with x ∈ . The concentration c(x, t + n ) of Equation (14) is the "initial condition" that has to be used for solving Equation (4) in the interval (t n , t n + 1 ). We remark that, when the vesicle concentration inside the synaptic regions is such that it is safely approximated by its mean value, i.e., N i /Vol( i ), for all i = 1, . . . , S, Equation (11) becomes Consequently, Equation (14) can be approximated as follows We compute the variation of the overall number of vesicles contained in the bouton by integrating Equation (14) over and using the properties of the characteristic functions (see Section 3.9), i.e., We remark that the number of vesicles does not vary, i.e., n − 1 is zero. In order to determine the solution to Equation (4) in each interval I m , with m = 1, . . . , N A + 1, we consider the restriction c m to I m , and discretize Equation (4) in time by using a first-order finite difference implicit Euler scheme. The result of this calculation reads where c m (x) ≡ c m x, t , with = {k, k + 1}, denotes the concentration in x ∈ at time t ∈ I m , and τ k = t k + 1 − t k > 0 is the amplitude of the kth time step. Then, we solve Equation (18) by using the Finite Volume Method. For a detailed presentation of the solution method we refer to standard literature on Finite Volumes. In this paper, we sketch very briefly the procedure to compute N i t − n numerically. To this end, we firstly cover the entire computational domain with a finite element triangulation. The nodes of the triangulation are denoted by x s , with s = 1, . . . , N n being the node index, and N n being the total number of nodes. The unknown of the problem, i.e., the vesicle concentration, is approximated by the sum c(x, t) ≈ N n s = 1 c s (t)ψ s (x), where c s (t) represents the time-varying value of concentration at the node x s , whereas ψ s is referred to as shape function. It is such that ψ s (x r ) = δ s r , with δ s r being the Kronecker-delta (it returns 1, if r = s, and it returns 0 otherwise). Secondly, we construct a dual grid by following the methods outlined in Frolkovic (1996). An "element" of the dual grid is referred to as finite volume. The finite volume associated with the node x s is denoted by B s in this work, and is obtained by joining the barycenters of the elements sharing the node x s with the midpoints of the edges connecting x s with its neighboring nodes. The finite volumes B s , with s = 1, . . . , N n , cover the entire computational domain. For details the reader is referred to Frolkovic (1996). The whole procedure is applied also to the synaptic regions i , with i = 1, . . . , S. Let then B ij specifically denote the finite volume associated with the node x ij ∈ i of the discretization covering the synaptic region i , with j = 0, . . . , M i being the nodal index of the discretization of i , and (M i + 1) the total number of nodes of this discretization. Then, it holds that Vol( i ) = M i j = 0 Vol(B ij ). Moreover, the number of vesicles N i (t − n ) is approximated as follows Finally, by substituting Equation (19) into Equation (9), we obtain the discretized form of the expression

NUMERICAL METHODS
The model equations were discretized with finite volumes in space and an implicit Euler method in time. The resulting systems of linear equations were solved with numerical multigrid techniques using Gauss-Seidel pre-and post-smoothers and LUdecomposition as base solver. The resulting fine grid problems consisted of approximately 20,000-40,000 grid nodes and were solved on parallel architectures, the HERMIT Supercomputer at the HLRS Stuttgart and on the Cekon and Cesari cluster of the G-CSC.

BUILDING THE COMPUTATIONAL DOMAIN
Bouton morphologies were constructed with Blender and ProMesh and discretized with a tetrahedral meshing algorithm using TetGen, that is incorporated in ProMesh (Reiter and Wittum, 2014). Examples of the 3D geometries can be found in Figure 1.

EMPIRICAL LAWS FOR THE IN SILICO RESULTS
From the simulations performed in Figures 5A,C,E we derived fitted empirical laws that reproduce the qualitative features of the measured post-synaptic response. The three parameters that were varied are bouton size, stimulation frequency and P o . In the following we denote the number of exocytosed vesicle as E(t, ω) and E(t, d) respectively, the stimulation frequency as ω, the vesicle output probability as P o , the number of synapses as S, the bouton diameter as d and use characteristic time points t i with i = 1, . . . , 6 for tuning the empirical laws. Note that the timepoints t 2 , t 4 , t 6 refer to the time-points in which the respective functions reach the value 0.01. The characteristic time points originate from heuristic fits using logarithmic representations resp. projections in x and y space allowing for heuristic determination of approximating curves. The same holds true for the E(t, . . . ) which are used in the following. All equations stated in the following are used in dimensionless form, yet for the sake of clarity we omit the procedure of making the equations dimensionless during the introduction of the following equations. In the following section we define θ (x) as

Bouton size-variation law
Considering the bouton size to be variable we can investigate the effect of the bouton size on the postsynaptic response that was simulated in Figure 5C. From the in silico experiments, we derive the following empirical law with S = 7, 10, 14 for d = 2, 3, 5 µm respectively and P o := 6%: with characteristic time points t 1 (d) = exp log (11.5) + log 465 11.5 · d 5 (23) t 2 (d) = exp log (80) + log 2900 80

Frequency-variation law
The effect of the stimulation frequency on the postsynaptic response was simulated in Figure 5E. From the in silico experiments, we derive the following empirical laws with S := 10 and P o := 6% with frequency-dependent characteristic time points t 4 (ω) = exp 8.8 − 0.775 log (ω)

P o -Variation law
The effect of the vesicle output probability P o on the post-synaptic response was simulated in Figure 5A. From the in silico experiments, we derive the following empirical laws with S := 10 and with the characteristic time points where

NORMS FOR QUANTIFYING THE SIMULATION RESULTS
In order to compare the experimental results with the simulation results, we applied two norms to the data: 1. Least-squares S, defined as S := 1 , with e i being the experimental data points and s i the simulated data points. p i denote the stimulus number.

SUMMARY OF THE MATHEMATICAL EXPRESSIONS USED IN THE MODEL
The symbols δ, and χ i denote, respectively, the Dirac's Delta distribution, the Heaviside function and the characteristic function associated with the set i . The symbol + is used to modify the definition of , as explained below. Given a set V ⊂ R d , with R d being the d-dimensional Euclidean space, the function χ V : R d → {0, 1} is said to be the characteristic function associated with V. It is defined as follows Let V ⊂ be a subset of and let f : → R be an integrable function over V. Then, the following identity holds true where f |V represents the restriction of f to V. In particular, if f is everywhere equal to unity, then it follows that The Heaviside function, , is defined as If the function is required to return unity also for ξ = 0, then the definition (34) is extended as follows The Dirac Delta is a distribution. Given an arbitrary function φ : R → R, the Dirac Delta satisfies the important property if φ is continuous at t ∈ R.

DISCUSSION
In this study we have developed a three-dimensional model of presynaptic vesicle dynamics to assess the physiological meaning of experimentally rather inaccessible bouton parameters like the bouton size with a given release site density, vesicle diffusion and the vesicle release probability. In combination with the stimulation frequency and an associated release of vesicles these parameters represent simplified but fundamental bouton output functions that should allow a systematic structure/function analysis of two very similar classes of glutamatergic boutons.

A THREE-DIMENSIONAL MODEL OF LOCAL VESICLE DYNAMICS
In addition to state of the art modeling approaches to synaptic or neuronal functions we used methods that represent the threedimensional morphology and biological processes derived from physical principles. Electrical signaling of neurons is often modeled as a one-dimensional process using the cable equation and compartmental one-dimensional space discretization (Hines and Carnevale, 2001). There are few approaches in cellular modeling that make use of 2D or 3D modeling techniques (Smart and McCammon, 1998;Ross et al., 2000;Meinrenken et al., 2002;Tai et al., 2003;Marpeau et al., 2009;Bielecki et al., 2010). Most commonly, stochastic models are employed to address the high complexity of three-dimensional simulations. A standard approach to in silico studies of presynaptic vesicle dynamics is to model the transition of vesicles between defined states with preconfigured transition-rate constants (Pan and Zucker, 2009), thus reducing the three-dimensional process of vesicle motion to a timedependent process (zero-dimensional in space). Considering the underlying question of our study presented here, we needed to develop a model that takes into account the spatial dynamics of vesicles in 3D. Modeling neurobiological systems requires incorporating the physical laws to which these systems abide. The use of deterministic models derived from continuum mechanics theory allows us to describe spatio-temporal vesicles dynamics based on physical laws which allows for a much more profound understanding of the problem compared to the use of heuristic models which are afterwards fitted to data. The parameters used for computations were taken either from the literature or EM sections and no parameter fitting was used. Combined with a three-dimensional representation of the morphology using unstructured grids, in this case different types of presynaptic boutons, this approach can answer questions regarding structure and function relationships. The equations are numerically approximated by a finite volumes approach and solved on parallel computers. The added mathematical techniques allowed us to investigate the interplay between structure and function at presynaptic terminals in a highly realistic manner. In addition to the three-dimensional representation of the underlying biophysics and three-dimensional morphology, our approach enables us to continually refine the model by including more of the systems properties, such as vesicle recycling, calcium dynamics and others depending on the current focus.

VESICLE DIFFUSION LIMITS THE AVAILABILITY OF PREEXISTING VESICLES AT ACTIVE ZONES
Our simulations modeled vesicle motility during stimulation as a simple undirected diffusion process. These simulations matched our experimental results (Figures 1E,F, 3, 4) suggesting that vesicle diffusion describes the motility of preexisting vesicles during periods of high frequency stimulation quite accurately. However, there is a debate whether or not active and directional transport mechanisms based on, e.g., the Myosin/Actin system are directly contributing to synaptic transmission (e.g., Ryan, 1999;Jordan et al., 2005;Tokuoka and Goda, 2006;Seabrooke et al., 2010;Kisiel et al., 2011;Seabrooke and Stewart, 2011). In particular in the Drosophila NMJ system recent evidence has implicated nonmuscle myosin II and the unconventional myosin VI in mobilizing synaptic vesicles during synaptic transmission (Seabrooke et al., 2010;Kisiel et al., 2011;Seabrooke and Stewart, 2011). In our study we focused on high frequency synaptic transmission in the absence of vesicle recycling and hence analyzed exclusively the motility of preexisting mature vesicles during extended periods of naturally occurring firing frequencies of 10-80 Hz. Under these conditions the modeled slow vesicle diffusion was sufficient to rather accurately mimic our eEJP-depletion experiments (Figures 3, 4) and they showed no obvious requirement for the incorporation of a faster vesicular transport mechanism. However, our results also showed that vesicle recycling plays a major role during prolonged high frequency stimulation as indicated by the much slower eEJP-depletion kinetics in the absence of BafA1 (Figure 3D). It seems therefore possible that fast active transport mechanisms are primarily involved in vesicle recycling processes that rapidly replenish active zones with vesicles during periods of prolonged stimulation. In contrast, preexisting vesicles that may for example be part of the reserve pool of vesicles and that have been assayed in this study seem to be recruited to active zones primarily by diffusion. Apparently, the glutamatergic boutons of Drosophila NMJs use both active and passive vesicular transport mechanisms. Whether and how they interact and what their relative contributions to the rather short natural activity patterns ( Figure 6C) are has to be analyzed in future studies. Based on our finding that vesicle diffusion is the primary factor limiting the replenishment of active zones with preexisting vesicles, we have assayed its effects on synaptic performance under various conditions. Depending on the activity of a given release site, which was determined by its release probability, its stimulation frequency and the typical number of stimuli per pattern, the local vesicle concentration dropped rapidly resulting in individual transmission failures and hence in a fast decline of the total bouton output (Figures 5A,B,E,F). The size of a bouton affects these parameters indirectly by offering, e.g., proportionally more active zones in larger boutons and simultaneously a larger functional reservoir of vesicles resulting in this example in a larger and more enduring bouton output, respectively (Figures 5C,D). Reliable transmission of a given activity pattern therefore depended strongly on the relative settings of these parameters. Based on these simulations we predict that a small type Is bouton with its smaller number of rather high P o -release sites is inappropriate for the transmission of long-lasting high stimulation frequencies ( Figure 5) and should quickly become diffusion-limited (white patches in right bouton of Figure 6B). This effect is governed by both bouton size and the corresponding number of release sites. In contrast, a large type Ib bouton with its larger number of low P o -release sites is likely to be better adapted to transmit such a pattern (Figure 5; left bouton in Figure 6B). Indeed, many experimental high frequency stimulation paradigms, in which both motor neurons were fired simultaneously by a common stimulation electrode showed a depression of early compound eEJP-amplitudes before reaching plateau values (Delgado et al., 2000).
Our predictions now revealed that the observed eEJPdepression in HFS-trains is likely due to an initial dominance and rapid decay of the Is output while the Ib output maintains its reliable transmission ( Figure 6A). It thus seems obvious that in vivo Is and Ib motor neurons can not be fired in a synchronous manner. Instead, it is more likely that the type-Is motor neuron is either fired with a low stimulus frequency or with only few stimuli, or both. Type Ib motor neurons are, in contrast, capable of reliably transmitting longer lasting high frequency patterns. We found that the natural firing activities of both motor neurons (Figures 6C,D) were strikingly consistent with the above predictions indicating that motor neuron firing patterns and NMJ-bouton properties are well-adapted features of this communication pathway.

THE PREVAILING FIRING PATTERN SHAPES BOUTON PROPERTIES
It has been reported in several central and peripheral systems including the Drosophila NMJs that long-lasting changes in the neuronal firing activity can change the physiological and morphological properties of the corresponding output synapses (Mollgaard et al., 1971;Hubel et al., 1977;Lnenicka and Atwood, 1989;Lnenicka et al., 2003;Sigrist et al., 2003). From these observations the concept emerges that the prevailing firing pattern of a neuron shapes the properties of its output synapses to ensure reliable transmission of its typical patterns. The simulations performed in this study revealed several tools a neuron can use to tune its output performance to its specific needs (Figure 7). Under these conditions the developing bouton has to find an economic and still effective/reliable compromise between its size, i.e., the number of vesicles and release sites and their release probabilities ( Figure 7B). A sustained high-frequency pattern would require a relatively large bouton with low release probability release sites. In contrast, if the prevailing activity pattern consists of only few APs at a rather low frequency an economic realization could be a smaller bouton with high-probability release sites. For a developing Drosophila NMJ these are the equivalents of type Ib and type Is boutons, respectively (Figures 1A, 7C). Early in larval development, when the innervated muscles are very small only few boutons of each type are required to ensure the transmission of the motor pattern and an efficient excitation-contraction coupling. The larger the muscle grows the less efficient excitationcontraction coupling becomes, requiring the proportional addition of more boutons of both types (Schuster et al., 1996). In addition to this developmental change of boutons, larval NMJs can undergo an experience-dependent bouton adaption (Sigrist et al., 2003). The morphology of a mature NMJ ( Figure 7C) therefore likely represents a cumulative readout of its prevailing firing pattern throughout development.

CONCLUDING REMARKS
Our three-dimensional functional model of multi-release site boutons allowed us to dissect the physiological meaning of experimentally inaccessible parameters like the bouton size, vesicle diffusion and the release probability. Based on the simulations with the three-dimensional model, we were able to derive an empirical and analytical description of the complex system subject to variations in bouton size, release probability and stimulation frequency. Our model sets the framework for a theoretical reconstruction and functional dissection of further realistic parameters such as vesicle recycling, the membrane potential and its dynamics during stimulation, the implementation of voltage-dependent calcium channels and intracellular calcium transients etc. Many of these parameters are experimentally difficult to access. Their implementation into a realistic theoretical model will therefore enable insights that have previously been impossible.