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

^{1}Bernstein Group Detailed Modeling of Signal Processing in Neurons, University of Heidelberg and University of Frankfurt, Heidelberg/Frankfurt, Germany^{2}Department of Simulation and Modeling, Goethe Center for Scientific Computing, University of Frankfurt, Frankfurt, Germany^{3}Bernstein Center for Computational Neuroscience Heidelberg-Mannheim, Heidelberg/Mannheim, Germany^{4}Department of Neurobiology, Interdisciplinary Center for Neurosciences, University of Heidelberg, Heidelberg, Germany^{5}Development Unit, European Molecular Biology Laboratory, Heidelberg, Germany^{6}Department of Mathematical Sciences, Polythecnic of Turin, Turin, Italy^{7}Department of Computational Neuroscience, Goethe Center for Scientific Computing, University of Frankfurt, Frankfurt, Germany

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.

## 1. 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-release-site 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 structure-function 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.

## 2. Results

### 2.1. 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., 2002, 2003) with an average active zone diameter of 300–400 nm (Meinertzhagen et al., 1998; Sigrist et al., 2002, 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.

**Figure 1. Functional reconstruction of synaptic boutons of Drosophila in silico. (A)** Confocal image of a

*Drosophila*larval neuromuscular junction (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).

**Figure 2. Two-dimensional sketch of the bouton model**. Region Ω_{0} denotes the bouton area in which vesicle motion is governed by diffusion. A cutout region in the center represents space occupied by organelles, such as mitochondria. Ω_{1} … Ω_{4} define synaptic active zones. In these subdomains, vesicles diffuse and are subject to exocytosis modeled by the term *f*_{syn}, that also includes the transition from the continuous to discrete scale in Ω_{1} … Ω_{4}. The organelle as well as bouton boundaries are treated as no-flux boundary conditions.

Synaptic active zones were represented by a multiedge-surface 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 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 a constant output probability *P*_{o} throughout a given simulation. For literature-based simulations, we used *P*_{o} = 29 and 7 % assigned to Is and Ib boutons respectively, values taken from Kurdyak et al. (1994); Cooper et al. (1995) instead of regulated *P*_{o} values that are known to be involved in synaptic facilitation and depression. This restriction enabled us to dissect the role of each of these parameters on presynaptic vesicle dynamics under otherwise constant conditions. Although vesicle recycling has an effect on the transmission efficacy, we did not incorporate vesicle replenishment into the model in order to avoid introducing still unknown variables. Since the experiments also focus exclusively on preexisting vesicles, model and experiments are comparable.

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).

### 2.2. 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 high-frequency 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 BafA1-treated 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.

**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.

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 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).

#### 2.2.1. 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 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.

**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.

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 stimulation-induced 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).

**Table 2. Comparing the experimental and simulated data at different frequencies and with the least-squares method S and L_{2}-norm (see Materials and Methods) shows good agreement between experimental and simulation data shown in Figure 4A**.

#### 2.2.2. 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 frequencies similar to the experimental results, we systematically varied each individual component using as standard settings *P*_{o} = 6%, *d* = 3 μm (10 active zones) and ω = 40 Hz, with initial vesicle density set to 400 μm^{−3}, in order to dissect their respective roles in determining the functional properties of different bouton types [Figure 5 (dotted)]. From these simulations we derived empirical laws [Figure 5 (lines)] as a heuristic fit of the *in silico* data (see Materials and Methods) in order to assist the interpretation of our simulations.

**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 %.

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%.

#### 2.2.3. 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).

**Figure 6. Rational design of a naturally performing NMJ in silico**.

**(A,B)**Simulation of the total

**(A)**or spatially resolved

**(B)**vesicle dynamics of a type Ib (3 μm, 10 active zones;

*P*

_{o}= 5%) and Is (2 μm, 7 active zones;

*P*

_{o}= 40%) bouton during continuous 40 Hz-stimulation. In both cases the initial vesicle density was set to 400 μm

^{−3}.

**(A)**Under these conditions the total output dynamics of combined Is and Ib boutons showed a characteristic biphasic depletion.

**(B)**Is boutons showed more depleted release sites than Ib boutons at similar stimulus numbers.

**(C)**Simultaneous recordings of the natural larval motor pattern in muscles 6 and 13 identified a synchronous firing pattern (arrows) originating from one motor neuron forming Is boutons on both muscles. The remaining asynchronous firing (arrowheads) originates from muscle specific motor neurons forming Ib boutons.

**(D)**Dissected frequency profile of natural larval motor patterns plotted as a function of the relative train length. Data represent means ± SEM from 10 identified Is- and 10 Ib-motor patterns from muscle 6.

**(E)**Simulation of the vesicle output of an

*in silico*NMJ (17 3 μm boutons,

*P*

_{o}= 7%; 19 2 μm boutons (7 active zones),

*P*

_{o}= 29%). 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)**.

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.

#### 2.2.4. 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, 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).

**Figure 7. Parameters defining the configuration of a NMJ**. **(A)** The prevailing firing pattern of a neuron sets the standard for the functional properties of its synaptic boutons. **(B)** In order to guarantee reliable transmission of a given firing pattern the already existing synaptic boutons may modify their genetically determined layout by adapting the here identified critical parameters such as the release probability and the bouton size (together with the number of release sites). **(C)** Throughout larval development individual boutons are added to the terminals of both motor neurons resulting in strings of boutons with typical morphologies.

## 3. Materials and Methods

### 3.1. Electrophysiological Recordings

Electrophysiological recordings of membrane potentials of body wall muscle 6 were performed on filet preparations of size-matched 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.

### 3.2. Parameters Used in the Model

#### 3.2.1. Bouton size

Paraformaldehyde fixed larval filet preparations were immuno-labeled 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.

#### 3.2.2. 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). (2) Ultrathin-sectioned larval boutons showed on average 1–2 vesicular profiles every 200 nm (Sigrist et al., 2000; Steinert et al., 2006) mounting up to about 500–1000 vesicles/μm^{3}. For our simulation experiments we have chosen a vesicle concentration between these estimates: 150–500 vesicles/μm^{3}.

#### 3.2.3. 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.

### 3.3. 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* ⊂ ℝ^{3} and an interval of time = [0, *T*] are considered. The total number of vesicles contained in *V* at time *t* ∈ 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 Ω ⊂ ℝ^{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:

Finally, we close the mathematical problem by prescribing the initial condition

where *t* = 0 denotes the origin of time.

### 3.4. 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}}{=}\frac{{n}}{{\omega}}$ ∊ , and

*n*= 1, …,

_{A}. In the following, we explain the terms used in the above equations:

1. Given the interval , that represents the lapse of time during which the system is observed, the natural number _{A} ∈ ℕ is the total number of times (modeled as points of ) in which the action potential stimulates the depletion of vesicles. The generic instant of time *t*_{n} ∈ , with *n* = 1, …, _{A}, corresponds to the *n*th stimulus, and is defined by ${{t}}_{{n}}{:}{=}\frac{{n}}{{\omega}}$, 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 *i*th 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 Θ(*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, φ : ↦ ℝ, 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*_{}_{A} + 1 = *T*, and write the interval = [0, *T*] as = ∪^{A + 1}_{m = 1}_{m}, where

In each open interval (*t*_{m − 1}, *t*_{m}), *m* = 1, …, _{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}}{=}\frac{{p}}{{\omega}}$ ∈ , with

*p*∈ {1, …,

_{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*(

**x**, t^{±}

_{p}) = lim

_{t → t±p}

*c*(

*), we obtain*

**x**, tThe “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}) d*v*, we obtain

We denote now by *c*_{n}, with *n* = 1, …, _{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, …,

_{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*_{Ω}(*t*^{+}_{n}) = *N*_{Ω}(*t*^{−}_{n}), if either Θ(*P*_{o} − *P*_{i}(*t*_{n})) or Θ_{+}(*N*_{Ωi}(*t*^{−}_{n}) − 1) is zero.

In order to determine the solution to Equation (4) in each interval _{m}, with *m* = 1, …, _{A} + 1, we consider the restriction *c*_{m} to _{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

**∈ Ω at time**

*x**t*

^{ℓ}∈

_{m}, and τ

^{k}=

*t*

^{k + 1}−

*t*

^{k}> 0 is the amplitude of the

*k*th 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}being the node index, and

_{n}being the total number of nodes. The unknown of the problem, i.e., the vesicle concentration, is approximated by the sum , 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}) = ⊔

^{Mi}

_{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 *h*(*x*_{s}, *t*^{−}_{n}):

### 3.5. 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 LU-decomposition 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.

### 3.6. 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.

### 3.7. 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 time-points *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

#### 3.7.1. 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

#### 3.7.2. 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

#### 3.7.3. *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 ω := 40 Hz

with the characteristic time points

where

### 3.8. 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}{:}{=}\frac{{1}}{{n}}\sqrt{{\displaystyle {{\sum}}_{{i}{\text{\hspace{0.05em}}}{=}{\text{\hspace{0.05em}}}{1}}^{{n}}{{d}}_{{i}}^{{2}}}}$

2. *L*_{2}-*norm*, defined as ${{L}}_{{2}}{:}{=}\sqrt{\frac{{{d}}_{{1}}^{{2}}\frac{{(}{{p}}_{{2}}{-}{{p}}_{{1}}{)}}{{2}}{+}{\displaystyle {\sum}_{{i}{\text{\hspace{0.05em}}}{=}{\text{\hspace{0.05em}}}{2}}^{{n}{\text{\hspace{0.05em}}}{-}{\text{\hspace{0.05em}}}{1}}{{d}}_{{i}}^{{2}}}\frac{{(}{{p}}_{{i}{\text{\hspace{0.05em}}}{+}{\text{\hspace{0.05em}}}{1}}{-}{{p}}_{{i}{-}{1}}{)}}{{2}}{+}{{d}}_{{n}}^{{2}}\frac{{(}{{p}}_{{n}}{-}{{p}}_{{n}{\text{\hspace{0.05em}}}{-}{\text{\hspace{0.05em}}}{1}}{)}}{{2}}}{{(}{{p}}_{{n}}{-}{{p}}_{{1}}{)}}}$

where ${{d}}_{{i}}{=}{\left|}\frac{{(}{{s}}_{{i}}{-}{\text{\hspace{0.17em}}}{{e}}_{{i}}{)}}{{(}{{s}}_{{i}}{+}{\text{\hspace{0.17em}}}{{e}}_{{i}}{)}}{\right|}$, with *e*_{i} being the experimental data points and *s*_{i} the simulated data points. *p*_{i} denote the stimulus number.

### 3.9. 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* ⊂ ℝ^{d}, with ℝ^{d} being the *d*-dimensional Euclidean space, the function χ_{V} : ℝ^{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* : Ω → ℝ 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 ϕ : ℝ → ℝ, the Dirac Delta satisfies the important property

if ϕ is continuous at *t* ∈ ℝ.

## 4. 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.

### 4.1. 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 three-dimensional 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 time-dependent 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.

### 4.2. 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 non-muscle 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 eEJP-depression 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.

### 4.3. 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 excitation-contraction 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.

### 4.4. 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.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgment

This work was funded by a BMBF grant Bernstein group DMSPiN (BMBF, 01GQ0803) to Christoph M. Schuster and Gabriel Wittum and a grant from the German Ministry of Education and Research (BMBF, 01GQ1003A) to Gillian Queisser and Christoph M. Schuster. Additional funding came from the NeFF Network of Frankfurt University to Gillian Queisser and Gabriel Wittum. We thank I. Heppner, M. Lampe, A. Nägel, S. Reiter and K. Xylouris, for helpful discussions. We would like to thank M. Stepniewski for help with the volume grid generation.

## Supplementary Material

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fncom.2014.00101/abstract

## References

Atwood, H., Govind, C., and Wu, C. (1993). Differential ultrastructure of synaptic terminals on ventral longitudinal abdominal muscles in drosophila larvae. *J. Neurobiol*. 24, 1008–1024. doi: 10.1002/neu.480240803

Bastian, P., Birken, K., Johannsen, K., Lang, S., Reichenberger, V., Wieners, C., et al. (2000). “Parallel solution of partial differential equations with adaptive multigrid methods on unstructured grids,” in *High Performance Computing in Science and Engineering*, eds W. Jaeger and E. Krause (New York, NY: Springer), 506–519.

Bielecki, A., Kalita, P., Lewandowski, M., and Siwek, B. (2010). Numerical simulation for a neurotransmitter transport model in the axon terminal of a presynaptic neuron. *Biol. Cybernet*. 102, 489–502. doi: 10.1007/s00422-010-0380-z

Bradacs, H., Cooper, R., Msghina, M., and Atwood, H. (1997). Differential physiology and morphology of phasic and tonic motor axons in a crayfish limb extensor muscle. *J. Exp. Biol*. 200(pt 4), 677–691.

Chouhan, A., Zhang, J., Zinsmaier, K., and Macleod, G. (2010). Presynaptic mitochondria in functionally different motor neurons exhibit similar affinities for ca^{2+} but exert little influence as ca^{2+} buffers at nerve firing rates *in situ*. *J. Neurosci*. 30, 1869–1881. doi: 10.1523/JNEUROSCI.4701-09.2010

Collins, C., and DiAntonio, A. (2007). Synaptic development: insights from drosophila. *Curr. Opin. Neurobiol*. 17, 35–42. doi: 10.1016/j.conb.2007.01.001

Cooper, R., Stewart, B., Wojtowicz, J., Wang, S., and Atwood, H. (1995). Quantal measurement and analysis methods compared for crayfish and drosophila neuromuscular junctions, and rat hippocampus. *J. Neurosci. Methods* 61, 67–78. doi: 10.1016/0165-0270(95)00024-O

Darcy, K., Staras, K., Collinson, L., and Goda, Y. (2006). Constitutive sharing of recycling synaptic vesicles between presynaptic boutons. *Nat. Neurosci*. 9, 315–321. doi: 10.1038/nn1640

Delgado, R., Maureira, C., Oliva, C., Kidokoro, Y., and Labarca, P. (2000). Size of vesicle pools, rates of mobilization, and recycling at neuromuscular synapses of a drosophila mutant, shibire. *Neuron* 28, 941–953. doi: 10.1016/S0896-6273(00)00165-3

Denker, A., Krohnert, K., and Rizzoli, S. (2009). Revisiting synaptic vesicle pool localization in the drosophila neuromuscular junction. *J. Physiol*. 587(pt 12), 2919–2926. doi: 10.1113/jphysiol.2009.170985

Dreyer, F., Peper, K., Akert, K., Sandri, C., and Moor, H. (1973). Ultrastructure of the active zone in the frog neuromuscular junction. *Brain Res*. 62, 373–380. doi: 10.1016/0006-8993(73)90699-9

Einstein, A. (1905). Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. *Annalen der Physik* 322, 549–560. doi: 10.1002/andp.19053220806

Fernandez-Alfonso, T., and Ryan, T. (2008). A heterogeneous resting pool of synaptic vesicles that is dynamically interchanged across boutons in mammalian cns synapses. *Brain Cell Biol*. 36, 87–100. doi: 10.1007/s11068-008-9030-y

Frolkovic (1996). “Finite volume discretization of density driven flows in porous media,” in *Finite Volumes for Complex Applications*, eds F. Benkhaldoun and R. Vilsmeier (Paris: Hermes), 433–440.

Gaffield, M., Rizzoli, S., and Betz, W. (2006). Mobility of synaptic vesicles in different pools in resting and stimulated frog motor nerve terminals. *Neuron* 51, 317–325. doi: 10.1016/j.neuron.2006.06.031

Gulyas, A. (1993). Hippocampal pyramidal cells excite inhibitory neurons through a single release site. *Nature* 366, 683–687. doi: 10.1038/366683a0

Hines, M., and Carnevale, N. (2001). Neuron: a tool for neuroscientists. *Neuroscientist* 7, 123–135. doi: 10.1177/107385840100700207

Hubel, D., Wiesel, T., and LeVay, S. (1977). Plasticity of ocular dominance columns in monkey striate cortex. *Philos. Trans. R. Soc. Lond. B Biol. Sci*. 278, 377–409. doi: 10.1098/rstb.1977.0050

Jan, L., and Jan, Y. (1976). L-glutamate as an excitatory transmitter at the drosophila larval neuromuscular junction. *J. Physiol*. 262, 215–236.

Jordan, R., Lemke, E., and Klingauf, J. (2005). Visualization of synaptic vesicle movement in intact synaptic boutons using fluorescence fluctuation spectroscopy. *Biophys. J*. 89, 2091–2102. doi: 10.1529/biophysj.105.061663

Kisiel, M., Majumdar, D., Campbell, S., and Stewart, B. (2011). Myosin vi contributes to synaptic transmission and development at the drosophila neuromuscular junction. *BMC Neurosci*. 12:65. doi: 10.1186/1471-2202-12-65

Kurdyak, P., Atwood, H., Stewart, B., and Wu, C. (1994). Differential physiology and morphology of motor axons to ventral longitudinal muscles in larval drosophila. *J. Comp. Neurol*. 350, 463–472. doi: 10.1002/cne.903500310

Kuromi, H., and Kidokoro, Y. (2000). Tetanic stimulation recruits vesicles from reserve pool via a camp-mediated process in drosophila synapses. *Neuron* 27, 133–143. doi: 10.1016/S0896-6273(00)00015-5

Lnenicka, G., and Atwood, H. (1989). Impulse activity of a crayfish motoneuron regulated its neuromuscular synaptic properties. *J. Neurophysiol*. 61, 91–96.

Lnenicka, G., and Keshishian, H. (2000). Identified motor terminals in drosophila larvae show distinct differences in morphology and physiology. *J. Neurobiol*. 42, 183–197. doi: 10.1002/(SICI)1097-4695(200005)43:2%3C186::AID-NEU8%3E3.3.CO;2-E

Lnenicka, G., Spencer, G., and Keshishian, H. (2003). Effect of reduced impulse activity on the development of identified motor terminals in drosophila larvae. *J. Neurobiol*. 54, 337–345. doi: 10.1002/neu.10133

Marpeau, F., Barua, A., and Josic, K. (2009). A finite volume method for stochastic integrate-and-fire models. *J. Comput. Neurosci*. 26, 445–457. doi: 10.1007/s10827-008-0121-7

Meinertzhagen, I., Govind, C., Stewart, B., Carter, J., and Atwood, H. (1998). Regulated spacing of synapses and presynaptic active zones at larval neuromuscular junctions in different genotypes of the flies drosophila and sarcophaga. *J. Comp. Neurol*. 393, 482–492. doi: 10.1002/(SICI)1096-9861(19980420)393:4<482::AID-CNE7>3.0.CO;2-X

Meinrenken, C., Borst, J., and Sakmann, B. (2002). Calcium secretion coupling at calyx of held governed by nonuniform channel-vesicle topography. *J. Neurosci*. 22, 1648–1667.

Mollgaard, K., Diamond, M., Bennett, E., Rosenzweig, M., and Lindner, B. (1971). Quantitative synaptic changes with differential experience in rat brain. *J. Neurosci*. 2, 113–127.

Pan, B., and Zucker, R. (2009). A general model of synaptic transmission and short-term plasticity. *Neuron* 62, 539–554. doi: 10.1016/j.neuron.2009.03.025

Reiter, S., and Wittum, G. (2014). Promesh – a flexible interactive meshing software for unstructured hybrid grids in 1, 2 and 3 dimensions. (in press).

Rollenhagen, A., and Lubke, J. (2006). The morphology of excitatory central synapses: from structure to function. *Cell Tissue Res*. 326, 221–237. doi: 10.1007/s00441-006-0288-z

Ross, M., Linton, S., and Parnas, B. (2000). Simulation studies of vestibular macular afferent-discharge patterns using a new, quasi-3-d finite volume method. *J. Comput. Neurosci*. 8, 5–18. doi: 10.1023/A:1008976030745

Ryan, T. (1999). Inhibitors of myosin light chain kinase block synaptic vesicle pool mobilization during action potential firing. *J. Neurosci*. 19, 1317–1323.

Satzler, K., Soehl, L., Bollmann, J., Borst, J., Frotscher, M., Sakmann, B., et al. (2002). Three-dimensional reconstruction of a calyx of held and its postsynaptic principal neuron in the medial nucleus of the trapezoid body. *J. Neurosci*. 22, 10567–10579.

Schuster, C. (2006). Glutamatergic synapses of drosophila neuromuscular junctions: a high-resolution model for the analysis of experience-dependent potentiation. *Cell Tissue Res*. 326, 287–299. doi: 10.1007/s00441-006-0290-5

Schuster, C., Davis, G., Fetter, R., and Goodman, C. (1996). Genetic dissection of structural and functional components of synaptic plasticity. i. fasciclin ii controls synaptic stabilization and growth. *Neuron* 17, 641–654. doi: 10.1016/S0896-6273(00)80197-X

Seabrooke, S., Qiu, X., and Stewart, B. (2010). Nonmuscle myosin ii helps regulate synaptic vesicle mobility at the drosophila neuromuscular junction. *BMC Neurosci*. 11:37. doi: 10.1186/1471-2202-11-37

Seabrooke, S., and Stewart, B. (2011). Synaptic transmission and plasticity are modulated by nonmuscle myosin ii at the neuromuscular junction of drosophila. *J. Neurophysiol*. 105, 1966–1976. doi: 10.1152/jn.00718.2010

Shtrahman, M., Yeung, C., Nauen, D., Bi, G., and Wu, X. (2005). Probing vesicle dynamics in single hippocampal synapses. *Biophys. J*. 89, 3615–3627. doi: 10.1529/biophysj.105.059295

Siddiqi, O., and Benzer, S. (1976). Neurophysiological defects in temperature-sensitive paralytic mutants of drosophila melanogaster. *Proc. Natl. Acad. Sci. U.S.A*. 73, 3253–3257. doi: 10.1073/pnas.73.9.3253

Sigrist, S., Reif, D., Thiel, P., Steinert, J., and Schuster, C. (2003). Experience-dependent strengthening of drosophila neuromuscular junctions. *J. Neurosci*. 23, 6546–6556.

Sigrist, S., Thiel, P., Reif, D., and Schuster, C. (2002). The postsynaptic glutamate receptor subunit dglur-iia mediates long-term plasticity in drosophila. *J. Neurosci*. 22, 7362–7372.

Sigrist, S., Thiel, P., Reiff, D., Lachance, P., Lasko, P., and Schuster, C. (2000). Postsynaptic translation affects the efficacy and morphology of neuromuscular junctions. *Nature* 405, 1062–1065. doi: 10.1038/35016598

Smart, J., and McCammon, J. (1998). Analysis of synaptic transmission in the neuromuscular junction using a continuum finite element model. *Biophys. J*. 75, 1679–1688. doi: 10.1016/S0006-3495(98)77610-6

Steinert, J., Kuromi, H., Hellwig, A., Knirr, M., Wyatt, A., Kidokoro, Y., et al. (2006). Experience-dependent formation and recruitment of large vesicles from reserve pool. *Neuron* 50, 723–733. doi: 10.1016/j.neuron.2006.04.025

Tai, K., Bond, S., MacMillan, H., Baker, N., Holst, M., and McCammon, J. (2003). Finite element simulations of acetylcholine diffusion in neuromuscular junctions. *Biophys. J*. 84, 2234–2241. doi: 10.1016/S0006-3495(03)75029-2

Tokuoka, H., and Goda, Y. (2006). Myosin light chain kinase is not a regulator of synaptic vesicle trafficking during repetitive exocytosis in cultured hippocampal neurons. *J. Neurosci*. 26, 11606–11614. doi: 10.1523/JNEUROSCI.3400-06.2006

Wagh, D., Rasse, T., Asan, E., Hofbauer, A., Schwenkert, I., Duerrbeck, H., et al. (2006). Bruchpilot, a protein with homology to elks/cast, is required for structural integrity and function of synaptic active zones in drosophila. *Neuron* 49, 833–844. doi: 10.1016/j.neuron.2006.02.008

Keywords: neuromuscular junction, boutons, modeling and simulation, structure-function relationships, morphology, firing pattern, adaption

Citation: Knodel MM, Geiger R, Ge L, Bucher D, Grillo A, Wittum G, Schuster CM and Queisser G (2014) Synaptic bouton properties are tuned to best fit the prevailing firing pattern. *Front. Comput. Neurosci*. **8**:101. doi: 10.3389/fncom.2014.00101

Received: 14 October 2013; Accepted: 07 August 2014;

Published online: 09 September 2014.

Edited by:

Misha Tsodyks, Weizmann Institute of Science, IsraelCopyright © 2014 Knodel, Geiger, Ge, Bucher, Grillo, Wittum, Schuster and Queisser. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Gillian Queisser, Department of Computational Neuroscience, Goethe Center for Scientific Computing, University of Frankfurt, Kettenhofweg 139, 60325 Frankfurt, Germany e-mail: gillian.queisser@gcsc.uni-frankfurt.de