- Research
- Open access
- Published:

# Stochastic models for the in silico simulation of synaptic processes

*BMC Bioinformatics*
**volume 9**, Article number: S7 (2008)

## Abstract

### Background

Research in life sciences is benefiting from a large availability of formal description techniques and analysis methodologies. These allow both the phenomena investigated to be precisely modeled and virtual experiments to be performed *in silico*. Such experiments may result in easier, faster, and satisfying approximations of their *in vitro/vivo* counterparts. A promising approach is represented by the study of biological phenomena as a collection of interactive entities through process calculi equipped with stochastic semantics. These exploit formal grounds developed in the theory of concurrency in computer science, account for the not continuous, nor discrete, nature of many phenomena, enjoy nice compositional properties and allow for simulations that have been demonstrated to be coherent with data in literature.

### Results

Motivated by the need to address some aspects of the functioning of neural synapses, we have developed one such model for synaptic processes in the *calyx of Held*, which is a glutamatergic synapse in the auditory pathway of the mammalia. We have developed such a stochastic model starting from existing kinetic models based on ODEs of some sub-components of the synapse, integrating other data from literature and making some assumptions about non-fully understood processes. Experiments have confirmed the coherence of our model with known biological data, also validating the assumptions made. Our model overcomes some limitations of the kinetic ones and, to our knowledge, represents the first model of synaptic processes based on process calculi. The compositionality of the approach has permitted us to independently focus on tuning the models of the pre- and post- synaptic traits, and then to naturally connect them, by dealing with “interface” issues. Furthermore, we have improved the expressiveness of the model, e.g. by embedding easy control of element concentration time courses. Sensitivity analysis over several parameters of the model has provided results that may help clarify the dynamics of synaptic transmission, while experiments with the model of the complete synapse seem worth explaining short-term plasticity mechanisms.

### Conclusions

Specific presynaptic and postsynaptic mechanisms can be further analysed under various conditions, for instance by studying the presynaptic behaviour under repeated activations. The level of details of the description can be refined, for instance by further specifying the neurotransmitter generation and release steps. Taking advantage of the compositionality of the approach, an enhanced model could then be composed with other neural models, designed within the same framework, in order to obtain a more detailed and comprehensive model. In the long term, we are interested, in particular, in addressing models of synaptic plasticity, i.e. activity dependent mechanisms, which are the bases of memory and learning processes.

More on the computer science side, we plan to follow some directions to improve the underlying computational model and the linguistic primitives it provides as suggested by the experiments carried out, e.g. by introducing a suitable notion of (spatial) locality.

## Background

During the last decades, the development of high-throughput technologies has produced a huge amount of information in the field of neurobiology, requiring the utilization of mathematical modeling to describe the complex dynamics of biological processes and stimulating new collaborations among biologists, physicists and computer scientists [1, 2]. The building blocks of neural systems are the neurons, which are specialized eukaryotic biological cells able to communicate with each other at highly specialized contact sites, called synapses. In general, each neuron consists of a somatic cellular body, on which a variable number of thin elongated structures, called dendrites, converge and from which a long single structure, called axon, emerges, branching in several synaptic terminals. The synaptic terminals of the transmitting neuron (the presynaptic element) send signals by releasing chemical molecules (neurotransmitters) to the dendritic, somatic or axonic, part of the receiving neuron (postsynaptic term) [3].

The synapses are the places of functional contacts between neurons, where the information is stored and transmitted from one to another neuron. Synaptic transmission is a complex process not completely understood: “Our current knowledge concerning synaptic transmission in neuronal networks of the brain is comparable to a puzzle in which most of the pieces are still missing” [4]. Current knowledge on synapses is based on the analysis of a limited number of experimental synaptic models, which were chosen for their experimental accessibility.

On the presynaptic side (sketched in Figure 1), one of the issues of the transmission concerns the involvement of calcium ions in the control of neurotransmitter release, which was hypothesized many years ago and more recently demonstrated [5]. The neurotransmitter release is mediated by exocitosis of synaptic vesicles (small elements containing the neurotransmitters) located at the presynaptic so-called active zones [6]. The electrical signals (action potentials) arriving at the synaptic terminal induce the opening of the *Ca*^{2+} channels. The transient elevation of the internal *Ca*^{2+} concentration in the presynaptic terminal triggers synaptic vesicle exocitosis, and hence the neurotransmitter release (calcium-triggered-release hypothesis). Interestingly, chemical messengers (intracellular) and modulators (extracellular) regulate the relationship between action potential and release in a synaptic terminal, which is also altered by the repeated activity. All these things make the presynaptic terminal a kind of computational unit, which changes its output based on its previous activity and ongoing modulation.

It must be noted that the *Ca*^{2+} signals were unaccessible to direct measurements up to few years ago. Theoretical and functional studies have suggested that calcium acts on presynaptic vesicles by a local, huge and short-lived elevation of its concentration. The locality and rapidity of the concentration variation render the study of this phenomenon not approachable with the conventional microscopic imaging techniques. Among the methods envisaged to overcome these limitations, one very fruitful approach is the so-called reverse approach, in which *Ca*^{2+} uncaging is induced in the presynaptic element. The uncaging method induces spatially homogeneous *Ca*^{2+} elevation, implying that measuring the *Ca*^{2+} fluorescent indicator gives an indication of the real *Ca*^{2+} sensed by the vesicles. This experimental method has been applied to the study of the large synapse of the auditory tract of the central nervous system, called *calyx of Held*. Moreover, it has been possible to build a minimal kinetic model for the process of the *Ca*^{2+} triggered vesicle release. Also, one can infer local *Ca*^{2+} signal waveform which is compatible with the experimental data on the time course and amplitude of release [7, 8].

Most of the models treating the issues of the calcium triggered release present some methodological limitations. These models, and among them the *calyx of Held* model, use differential equations to describe the time course of [*Ca*^{2+}], the *Ca*^{2+} concentration (*mole* × *liter*^{−1}) interacting with the synaptic vesicles. This approach implies that [*Ca*^{2+}] is continuous, while it is clearly not [9]. For example, a cubic volume with edge length 60 *nm* and a *Ca*^{2+} concentration of 10 *μ* M contains about a single free ion. Another common assumption is that the binding of the *Ca*^{2+} to the release sensors of the vesicle does not affect the [*Ca*^{2+}] concentration [9]. Also this assumption is not exactly adequate: considering that the dimensions of the vesicle diameters range in the interval 17-45 *nm*, an average interaction volume, a cube with edge length 60 *nm* say, would contain few *Ca*^{2+} ions, and when some of them bind to the vesicle sensors, the number of calcium ions could change substantially.

In general, when the fluctuations in the molecular population levels are relevant, e.g. when the numbers per unit volume of the molecular species involved are small, the stochastic approach has to be preferred. The fluctuations, that can be seen as random variations about a mean number of molecules, can significantly alter the dynamics of biochemical pathways. Just to cite few examples, fluctuations might decrease or increase the steepness of a nonlinear stimulus-response relationship. In other cases, via the so-called stochastic resonance, they can increase the reliability of response to small signals [10].

Hence, in these cases, the use of a stochastic approach appears to be much more appropriate, since the deterministic approach fails to capture the nature of chemical kinetics, which at low concentrations is discrete and stochastic [11]. In the stochastic approach, the system is described by the so-called “master equation”, which is usually intractable. A stochastic simulation algorithm has originally been proposed in [12] to overcome this difficulty.

On the postsynaptic side (sketched in Figure 1), the signal is transmitted when the neurotransmitters and the membrane receptors interact. The synaptic function is determined by many factors, such as the the speed of clearance of the neurotransmitters from the cleft, which determines the amount of postsynaptic receptors that bind the transmitters. This mechanism is also still poorly understood, because these events are beyond the power of experimental observations [13]. Due to its considerable size, the *calyx of Held* synapse represents a useful model system for studying the postsynaptic element, too. For this reason it has been the object of several experimental and theoretical works, e.g. [14, 15]. In general, the description of the postsynaptic interaction is made through deterministic approaches, which mainly rely on the use of molar concentration for the quantitative measure of the neurotransmitters [14, 16]. However, these suffer from the same limitations discussed for the case of the presynaptic vesicle release, and make the use of stochastic models preferable.

It is worth noting that the approach to devise models of the whole synaptic function has been mainly based so far on deterministic ODE descriptions [16] and the utilization of a stochastic approach, typically based on [12], is just at the beginning [17].

Recently, stochastic techniques have also been adopted in computer science to model quantitative aspects of interactive systems within concurrency theory. The aim of concurrency theory is to model the behaviour and the structure of systems composed of autonomous computational entities, which dynamically interact with one another, possibly reconfiguring the system itself. At the beginning, stochastic models have been used to study performance/time related properties, e.g. [18]. The strong analogies between concurrent and living systems, summarised in “cells as computation” [19], are contributing to the development of Systems Biology [20, 21], a systemic approach to living system modeling. According to this metaphor, cells, molecules and biological “active” components, i.e. those capable of exhibiting a behaviour, are assimilated to computer processes, the computational units of a concurrent software system. By communicating, processes may exchange information or synchronise themselves. Then, biological interaction corresponds to process communication. A biological experiment, or biological activity in general, has then a direct correspondence into a computation. That is, biological processes can not only be simulated by *in silico* experiments, but also it is possible to formally reason about their computational models and infer properties of interest. *Process calculi* are a formalism to describe such models: systems are *compositionally* described in terms of suitable *abstractions* of their component behaviour. Several process calculi whose “operators” are oriented to describing different aspects of biological interaction have been proposed, e.g. [22–25]. Some of these calculi have been equipped with a stochastic semantics in order to study the quantitative dynamics of systems, e.g. [23, 26–28]. This approach benefits from conjugating the abstract and compositional algebraic models, the possibility of precisely describing their semantics and formally reasoning about them, and the quantitative analysis provided by stochastic semantics. Executable implementations of the calculi and analysis tools are provided.

In this context, motivated by addressing some aspects of the functioning of neural synapses, we have first developed a stochastic model of the calcium triggered release in the *calyx of Held* synapse. Our work starts from a deterministic model presented in [7], from which we have derived a suitable stochastic model. This has subsequently been formalised in terms of Pi-calculus [29], a process calculus designed for concurrency and mobility. The behaviour of *Ca*^{2+} and vesicles has been described and composed to form the presynaptic terminal. Then, by the same approach, we have developed a stochastic model of the postsynaptic receptor kinetics, starting from the semi-deterministic model presented in [14, 16].

Finally, we have obtained a whole model of signal transmission throughout the synapse by suitably composing the two models of the pre- and post- synaptic sides developed individually. It must be noted that our model represents a fragment of the whole synapse focusing on a cluster of functionally related active zones. However, due to the morpho-functional organization of this synapse consisting of similar, independent active zones functioning in parallel, the behavior of the whole synapse can be easily obtained by scaling operation.

Model development has benefited from the above mentioned features, like modular design, abstract representation of the component functioning and stochastic interpretation of the system dynamics. Modular and abstract design, in particular, has been revealed to be suitable to facilitate the incremental development of a more complete and detailed model.

*In silico* experiments have been carried out by means of the stochastic Pi-calculus simulator SPiM [30], which represents one of the most complete and expressive simulation environments for stochastic calculi currently available. Obtained results are coherent with those from the deterministic models and with others in literature. The expressiveness of the model has allowed us to easily describe environmental conditions of the experiments, such as modulating presynaptic calcium waves. This has helped us to identify some critical points of the model and the parts of the overall process that seem to play a more relevant role. Beyond the experiments here reported, our approach appears to be interesting in the long term for the development of more comprehensive stochastic models of synaptic functioning. Moreover, the experiments carried out so far have been given interesting insights about the linguistic mechanisms that could usefully extend the expressiveness of the calculus.

The model is described in the next section, where we also report on some of the experiments done. The computational implementation of the model is reported in the last section *Methods*. Perspectives of the approach are discussed in *Conclusions*.

## Results and Discussion

We have applied a stochastic approach, based on the algorithm introduced in [12], to describe the presynaptic calcium triggered release mechanisms and the postsynaptic channel opening, studied in the model system of the synapse *calyx of Held*. For the presynaptic terminal we have modeled both the cases of step-like and wave-like calcium uncaging, while we have considered membrane receptor activity for the postsynaptic side.

### Step-like calcium uncaging

Our starting point was the following phenomenological kinetic model, described in [7]. (Slightly imprecisely, hereafter we shall use the notation A+B\begin{array}{c}\stackrel{to}{\to}\\ \underset{back}{\leftarrow}\end{array}C+D\begin{array}{c}\stackrel{\dots}{\to}\\ \underset{\dots}{\leftarrow}\end{array}\dots to represent the reversible reaction from *A* and *B* to *C* and back, and then, the produced *C* is involved, together with *D*, in the next (reversible) reaction producing …. I.e., the formula does not intend to represent neither *A* + *B* ← *C* + *D* nor *A* + *B* → *C* + *D*). Vesicles are activated through a process of five calcium binding steps and eventually released with a given kinetic rate constant. This process is determined by a cooperativity factor *b*. Vesicles are represented as *V*, the intracellular calcium as C{a}_{i}^{2+}, the released vesicles as *T* and the kinetic rate constant relative to vesicle release is γ = 6000 *s*^{−1}. The values of the other constants of the model are *k*_{
on
} = 9 × 10^{7}*M*^{−l}*s*^{−l}, *k*_{
off
} = 9500 *s*^{−1} and *b* = 0.25:

The parameters of the kinetic model were computed by a fitting of the experimental data. These had been obtained by elevating the intracellular presynaptic [*Ca*^{2+}] in a controlled, homegeneous and step-like manner [7].

The previous equations have been transformed by utilising the relationship between the stochastic rate constants (*c*) and deterministic rate constants (*k*) [31]. For reactions of the first order, *c* = *k*. For reactions of second order, the relationship becomes *c* = *k*/(*NA* × *Vol*), where *NA* represents the Avogadro's number and *Vol* the volume of the reaction. Hence, in order to determine the values of the stochastic rate constants, we need to estimate the value of *Vol*. Spatially, the calyx of Held is organized as a “parallel” arrangement of a large array of active zones, ranging from 300 to almost 700 [32]. Each active zone contains up to 10 vesicles and they are clustered in groups of about 10 of them, in a volume having a diameter of almost 1 μm. Each action potential activates all the active zones. Such particular morpho-functional organization of this synapse has allowed us to model a subunit of the presynaptic element, consisting of a cluster of 10 active zone, each containing 10 vesicles, in a volume of 0.5 10^{−15}*liter*. With this volume estimate, we have obtained the following values for the stochastic constants: *b* = 0.25, *c*_{
on
} = 9 × 10^{7} / (6.02 × 10^{23} × 0.5 × 10^{−15}) *s*^{−l} = 0.3 *s*^{−l}, *c*_{
off
} = 9500 *s*^{−l}, γ = 6000 *s*^{−l}, and the following numbers of *Ca*^{2+} ions: 300, 3000 and 6000, corresponding to molar concentrations [*Ca*^{2+}] of 1, 10 and 20 *μ* M. The equations of the stochastic model follow:

Note that, according to [7], we represent the fusion of each vesicle as the production of one released vesicle *T*, abstracting, for simplicity, from the actual quantities of neurotransmitters involved in the final step of the release process.

We have performed a series of simulations which confirmed the results obtained with the deterministic model: high sensitivity of vesicles to calcium concentrations [7]. Moreover, it is known that the local calcium concentration can be much lower than 100 *μ* M in the calyx of Held [7], while several other synapses require a calcium concentration in the range of 100-300 *μ* M for triggering vesicle releases [6]. Our results also confirm this fact, by showing that concentrations as low as 1, 10 and 20 *μ* M are able to deplete the releasable pool in a few milliseconds. In the following we will report a sample of our results. Each one consists of three pictures reporting the time course of *Ca*^{2+}, the logarithmic scale of it, which shows intermediate states of calcium binding, and a focus on vesicle activation (**Vstar**) and release (T), respectively. Figure 2 shows simulation results for the following parameters: *V* = 100; *Ca*^{2+} = 6000; *c*_{
on
} = 0.3; *c*_{
off
} = 9500; γ = 6000; *b* = 0.25. It can be observed that the pool of vesicles is 80% depleted within 3 *ms*, coherently to the experimental findings [7].

### Wave-like calcium uncaging

The technique of step-like calcium uncaging overcomes the limitations of traditional imaging methods [7] in resolving the local (and rapidly decaying) *Ca*^{2+} signal triggering the vesicles release. The step-like *Ca*^{2+} uncaging has shown a high sensitivity of the synapse under consideration to a uniform and sustained *Ca*^{2+} elevation in the range of 10 *μ* M (while, otherwise, *Ca*^{2+} concentration must be in the range of 100-300 *μ* M to have an effect). Yet, we needed to produce a spatially uniform rapidly decaying *Ca*^{2+} transient (wave-like calcium uncaging) in order to ascertaining whether very short *Ca*^{2+} elevations, as those following the action potential invasion of the synaptic terminal, suffice to induce vesicle release. Hence, the step-like uncaging has been modeled in order to measure the sensitivity of vesicles to uniform and sustained calcium concentration, whereas the wave-like uncaging has been modeled in order to ascertain the sensitivity of vesicles to short-lived variations of calcium concentration.

The experiments and models on *Ca*^{2+} uncaging [7] showed a high sensitivity of vesicle release in response to a uniform elevation of [*Ca*^{2+}] in the range 10 *μ* M. It was unclear whether very short [*Ca*^{2+}] elevations suffice to induce a release similar to that elicited during an action potential. A recent experimental work [33] has addressed this issue. A spatially uniform and very rapidly decaying [*Ca*^{2+}] transient was induced in the presynaptic element of a calyx of Held synapse by *Ca*^{2+} uncaging in the presence of added *Ca*^{2+} buffers. Remarkably, the short-lived elevation of calcium concentration triggers vesicle release. We have introduced in our model a simple mechanism of calcium extrusion utilized in a previously developed model [16], adapting the rate constants to fulfill our needs:

where C{a}_{o}^{2+} is the extruded calcium, *P* is an abstraction of a pumping mechanism, *c*_{1} = 8 *s*^{−1}, *c*_{2} = 25 *s*^{−1} and *c*_{3} = 10000 *s*^{−1}. Starting with 1000 *P*, we have obtained a simulated calcium wave lasting about 1 *ms* and with a half width 0.5 *ms*, conforming to the experimental requirements [33]. The left side of Figure 3 shows a calcium wave lasting about 1 *ms* with a half width of about 0.5 *ms* and a peak value of about 6000, corresponding to a peak calcium concentration of 20 *μ* M. The right side of the same figure depicts the release of one vesicle. Considering that a whole presynaptic element can be made of about 70 of our simulated clusters, this implies that a single action potential, and accordingly a single calcium wave, is able to release a significant amount of vesicles. This is also along the line of the experimental findings [7, 8, 33].

### Membrane receptor activity

We have then focused on the postsynaptic side of calyx of Held synapse and we have modeled the part of the postsynaptic membrane functionally related to the presynaptic mechanisms previously described. Our model considers postsynaptic membranes consisting of AMPA receptors and has been derived from a previously developed deterministic model, see [14, 34] for details.

The equations of our stochastic model are:

Slightly abusing the notation, the symbol *T* represents here the neurotransmitter molecules that bind to the receptor-gated channels, while in the presynaptic model it has been used to represent the whole content released by a single vesicle. Receptor-gated channels *C*_{0} are activated through a process of two reversible *T* binding steps: from *C*_{0} to *C*_{1}, intermediate channel states, and then to *C*_{2}, the activated channels. Either the activated channels evolve to the states *O*_{1} and *O*_{2}, open to the ion flux transmitting the signal (with different kinetics); or the activated channels evolve to the desensitized states *D*, which are not permeable to the ion flux. In the second case, the fraction of channels that open during a synaptic response decreases. All transitions are reversible.

The synaptic response, on the postsynaptic side, is determined by the time course of the rising phase of the synaptic ion current. This is a function of the opening rate of the receptor-gated channels and of the neurotransmitter concentration. We monitor the numbers of open *O*_{1}, *O*_{2} and densensitized *D* channels.

It is worth noting that the number of the *T* stochastically assumed to be released in the synaptic cleft has been determined coherently with the amount stochastically produced by a single calcium wave in the presynaptic terminal in our experiments (about 3 to 6 vesicles). Each “presynaptic” *T* corresponds to a wave of neurotransmitter whose duration and amplitude has been calculated on the basis of several hypotheses, below discussed.

The values of the kinetic rates are: *r*_{
b
} = 400 *s*^{−1}, {r}_{u}^{1} = 6 *s*^{−1}, {r}_{u}^{2} = 86000 *s*^{−1}, *r*_{
d
} = 900 *s*^{−1}, *r*_{
r
} = 64 *s*^{−1}, {r}_{o}^{1} = 100000 *s*^{−1}, {r}_{o}^{2} = 2000 *s*^{−1}, {r}_{c}^{1} = 2000 *s*^{−1}, {r}_{c}^{2} = 250 *s*^{−1}. These figures have been determined from the deterministic models in [5, 14], accordingly to the transformation already described in the section Step-like calcium uncaging.

In the definition of *r*_{
b
}, the “activation” rate of the receptor, a further difficulty arose: in order to determine the values of the stochastic rate constant, we not only need to estimate the *volume* occupied by the neurotransmitter in the cleft, but also the effective *interaction volume* where the interactions happen against the postsynaptic *surface*. We have computed the *volume* occupied by the neurotransmitters by knowing the number of them contained into a single vesicle (6000-8000) and the value of their concentration when they have spread into the synaptic cleft: 10^{−3}*M*. We obtained *volume* = 10^{−17}*liter*. Following [13], we estimated the effective *interaction volume* by dividing *volume* by 1/160, obtaining the result of 6 10^{−20}*liter*. We adopted this value for the conversion of the deterministic rate constant into the stochastic one, obtaining *r*_{
b
} = 400. Finally, the *interaction volume* determined permits us to estimate the number of neurotransmitters involved: about 50 *T*, as peak value. The value of *C*_{0}, which represents the number of membrane channels, has been assumed to be 100 according to known data, see e.g. [5].

We have performed a series of simulations and we report a sample of our results in Figure 4, which consists of three pictures showing the time course of the neurotransmitter *T* and the channel states *O*_{1}, *O*_{2}, *D*, *C*_{1}, *C*_{2}, *C*_{0}; the logarithmic scale of them; and a focus on the time course of *O*_{1}, *O*_{2}, *D*. It can be noted the very short duration of the *T* waves and a slower kinetics of the postsynaptic activated channels (compared to the duration of the *T* waves). This corresponds to the prolonged time course of the synaptic currents experimentally observed [34, 35]. Moreover, it can be observed that for each *T* wave, i.e. for each vesicle released according to our assumptions, the number of open channels is of the order of tens, fitting with experimental observations [5]. Finally, the buildup of the desensitization is also evident, again coherently with experimental data.

### Complete synapse signal traversal

A comprehensive model of the whole synapse has been devised by building upon the independently developed models of the pre- and post- synaptic processes, which have been considered in isolation in the previous sections. The calcium pulse induces the release of vesicles, represented abstractly in the presynaptic model as the release of a *single* neurotransmitter *T* by each vesicle. In the postsynaptic model, each *T* has been interpreted as a wave of neurotransmitters, with amplitude and duration determined as above. It must be underlined that in the study of synaptic transmission many of the elements we considered are generally neglected, such as the stochasticity of the vesicle release and the calculation of the real amount of neurotransmitter sensed by the postsynaptic receptors.

In Figure 5, we report a sample of results about the simulation of the model regarding a calcium wave with the same mechanisms and coefficients of those in the section about wave-like uncaging. The left part displays the calcium wave, amongst other values. The logarithmic scale of the same picture, in the middle, highlights the dynamics: two vesicles are released within 2ms and the corresponding two waves of neurotransmitters are clearly visible (also clearly distinguishable in the right part, the *T* curve). These waves start processes that lead to the opening of some receptor channels in a very short time and lasting for few milliseconds longer. This is coherent with experimental data, as described in the section about membrane receptors. It is worth noting that the *T* waves we supplied when describing the membrane receptor dynamic in isolation, are now induced stochastically by the calcium wave.

Beyond validating the consistency of the model, some preliminary experiments have also been run about open issues in short-term plasticity (reported in the subsequent section about plasticity). Indeed, the interaction between presynaptic and postsynaptic mechanisms generates short-term plasticity events in a way not yet fully understood [15].

### Parameter sensitivity analysis

We performed a parameter variation study (sensitivity analysis) in the case of step-like and wave-like calcium uncaging and for the postsynaptic terminal. The sensitivity analysis is a standard procedure in modeling and is performed with the aim of evaluating the robustness of the model, i.e. the tolerance of the release process to parameter variation, and to identify, if any, critical parameters.

For the presynaptic terminal, we have run simulations for different values of the number of vesicles (reference value 100, other values: 10, 50, 200 and 500), the number of calcium ions *Ca*^{2+} (reference values 300, 3000 and 6000, other values: 12000, 18000, 24000), the stochastic coefficients *c*_{
on
} (reference value 0.3, other values: 0.1, 0.2, 0.4 and 0.5), *c*_{
coff
} (reference value 9500, other values: 5500, 7500, 11500 and 13500), *b* (reference value 0.25, other values: 0.1, 0.2, 0.3 and 0.4) and γ (reference value 6000, other values: 2000, 4000, 8000 and 10000). One of the results of this analysis is that the forward coefficient *c*_{
on
} seems to have a critical role: higher values of this coefficient correspond to a faster release, in the step-like case, and to a switch from no-release to a consistent release, in the wave-like case. A selection of analysis results is reported below.

For the step-like case, Figures 6 and 7 show runs for *c*_{
on
} = 0.5 and *b* = 0.4, respectively. In Figure 6 we observe a high increase of the release rate, while in Figure 7 we observe a lower and more uniform release rate, both compared to the reference case (Figure 2). A simple variation of these coefficients changes the dynamics of the release processes in an unpredictable manner, producing in one case a variable and decreasing rate of release (see the change in the slope of the curve in between 1 and 1.5 *ms* in Figure 6) and in the other case a stationary rate of release (Figure 7). This kind of experiments are of interest when addressing the problem of the variations in the release rate, one of the still obscure phenomena which have been observed about vesicle release. These variations have been explained by the recruitment of new vesicles within the same active zone or by a different sensitivity to calcium ions of the vesicles belonging to the same cluster. Further investigations with our approach are ongoing about these controversial experimental data [7, 8, 33].

For the wave-like case, Figures 8 and 9 show runs for *c*_{
on
} = 0.1 and *c*_{
on
} = 0.5, respectively. It is interesting that in Figure 8 no release appears, whereas in Figure 9 a significant release (up to 8 vesicle in 1 *ms*) can be observed. These two different behaviour are obtained with the same (shape, peak value and duration) calcium wave. Hence the same very short lived calcium wave can be not enough or, on the contrary, even too much for inducing vesicle release. This is certainly a significant finding for the interpretation of recent experimental results, in which most of the attention has been directed, so far, only to the issue of calcium concentration variation [33]. A modulation of this coefficient might be a way for the synapse to store information, increasing or reducing the vesicle release rates (and ultimately it might be relevant for memory and learning).

We have also performed sensitivity analysis of parameters for the postsynaptic model. While the simulated channel dynamics appears overall suitably coherent with the experimental findings, the understanding of the parameter variation results is still under investigation. Amongst the other performed experiments, we report on the variation of two parameters. These two experiments are sufficient in order to illustrate the great variability of the roles of the postsynaptic parameters, which, indeed, are still poorly understood.

Figure 10 illustrates the simulation results obtained by varying of the coefficient {r}_{u}^{1} from 6 to 60. This is relevant for the unbinding rate of the transmitter-channel complex. This variation of one order of magnitude does not seem to significantly alter the kinetics of the main opening channel mechanism, represented by *O*_{1} (compare the rightmost pictures of Figure 10 and Figure 4). The slight variations in the kinetics of *O*_{2} and *D* could be explained in terms of stochastic fluctuations.

Figure 11 illustrates the simulation results obtained by varying the coefficient *r*_{
d
} from 900 to 9000. This represents the rate of desensitization. In this case the effects of a variation of the same order of magnitude as before are more marked: the time course of *D* is clearly faster and this impacts on the peak capability of *O*_{1}, too.

### Plastic events throughout the complete synapse

Some experiments addressing plasticity events in the whole synapse have been carried out using the complete synapse model. Besides describing some events of short-term plasticity, our model helps distinguishing the presynaptic and postsynaptic influence in these processes. It is know that in the process of short-term synaptic depression, distinct roles are played by vesicle depletion and postsynaptic receptor desensitization. It has also been shown that desensitization has a role in the synaptic depression only when the frequency of the action potentials is above 10 Hz [15]. Figure 12 shows results about an experiment regarding a train of calcium waves at a frequency of about 100 Hz, which mimics action potentials like neural signals. The left part displays calcium waves, and other values, shown in the central part in logarithmic scale. This makes clear the dynamics of the many processes involved, e.g. that of the neurotransmitter and the number of open channels. The right part shows a selection of the most significant values. Remarkably, a buildup of desensitization is visible starting from the fourth wave (see the increase of value *D* in between the third and fourth wave). Note, also, that the same amount of neurotransmitter causes different effects as far as the number of opened channels is concerned (compare the ratio between *O* 1 and *T* in the first two waves with the same ratio in the sixth and the seventh wave). One advantage of our approach is the possibility of quantifying the number of vesicles (recall that each vesicle correspond to a *T* wave with a peak of about 40 neurotransmitters and it is possible to see superpositions of some of them) stochastically associated to each action potential and of measuring at the same time the number of receptor channels open or desensitized. This means that it is possible to distinguish effects due to vesicle depletion, when it builds up, from the effects due to receptor channels.

## Conclusions

Our results are encouraging about the validity of the stochastic approach in studying the synaptic process, which consists of many discrete-like events and involves, for instance, arrays of vesicles and many different molecules. In the process of synaptic transmission, many of these molecules play roles still not fully understood [6]. We have studied parts of the process by adapting data from a number of experimental models. In doing this, we have overcome some limitations of several of these models, e.g. the assumption that concentrations are controlled and homogeneous, so that, for instance, the issue of spatial locality can be not considered. Our discrete approach has instead allowed us to address locality-related issues, like determining the actual interaction surface in the synaptic cleft, according to the quantity of released neurotransmitters and the volume they occupy. This approach appears amenable to also describe other “spatial” aspects of interest.

We built our model in a compositional manner. We first developed a presynaptic model and a postsynaptic model. Once tuned and validated against available experimental data, these models have been integrated to describe the signal traversal throughout the whole synapse. The compositional properties of the chosen approach play a crucial role here: designing and implementing a thin interface was sufficient to obtain the executable model for the whole synapse, with no change in those for the pre- and the post-synaptic processes.

Our model can be further refined along several directions. For instance, the use of a single uniform cluster of releasable vesicles is correct when the process is studied in a short interval of time. When events take place during long-lasting neural activity, more details on vesicle trafficking and cluster compartimentalisation should be introduced. Moreover, the vesicle release is regulated by many intracellular signal pathways, which influence the number and the speed of recruitment of the release-competent vesicles [36].

Embedding these aspects in the model might shed further light on the ways the nervous system processes and stores information. In order to support these, and others, developments the underlying process calculus needs extensions, too. Surely, the problem of expressing locality, already tackled by several calculi, could be valuably addressed within a stochastic viewpoint, together with the mentioned issue of model interfaces.

## Methods

In order to provide a grasp of the synapse formal model based on the Pi-calculus, we briefly sketch some of its parts regarding the presynaptic trait. Some excerpts from our model are reported in Figure 13. The language models interaction as pairs of input/output actions over the same communication channel (?c/!c). These atomic actions can be composed in a sequence (!c;?d; !c) or in alternative choices (?c or ?d) so as to form a possibly parametric process (p(…)= …). Processes can run in parallel (p() |q()). Initially, a command sets the duration of the simulation (here 0.005*s*), then some of the stochastic parameters are defined. Communication channels can be (dynamically) created by the new command and have associated a stochastic rate (this and the current quantities of reactants determine the probability of a reaction to “happen” through the channel). In our model one calcium ion, represented by the process ca() without any parameter, can interact with a vesicle v() over the channel vca with rate con5=1.5 (beyond being able to do other things). After this communication, ca() disappears (it becomes the null process ()) and v() becomes v_ca(), which represents the binding of the two. This realises a second order reaction. First order reactions are modeled as interactions with a *single* dummy molecule (so as not to alter stochastic dynamics). For instance, v_ca() can then either accept other calcium bindings, engaging in a second order reaction, or degradate back to an unbound vesicle v() by communicating with the dummy Dv_ca() through bvca, in a first order reaction. Dv_ca() restores ca() and itself.

So far, one can easily realise how the system has been described by specifying simple atomic behaviour, basically corresponding to chemical reactions, and then composing them together. The parametric process w(cnt:int) allows us to suitably modulate the calcium wave. After a stochastic delay, if its integer parameter cnt is positive then it replicates 80 copies of itself, with the parameter decreased, and 80 parallel copies of ca(), otherwise it dies. This realises an exponential growth, which can be controlled in its rapidity and quantity by the delay rate and the parameter. Finally, the initial state can be populated specifying how many molecules of each specie are present (one calcium wave, 1000 pumping molecules, 100 vesicles and the needed dummy molecules). The SPiM simulator produces the dynamics shown in the Figures 2,3,4,5,6,7,8,9,10,11,12.

It is also worth underlining the advantages provided by the compositionality of the approach, well shown by the construction of the complete synapse model. There, two models about different processes (happening in separate volumes), independently developed and tuned, have been connected by suitably mapping boundary processes, i.e. by defining a suitable interface between the models of the pre- and the post-synaptic parts. What was abstracted as the *T* single neurotransmitter released by a vesicle has become in the complete model a suitable wave of release. Actually, the presynaptic t() process now becomes the tt(1) process, which is the wave generator used in the postsynaptic model to spawn the expected amount of neurotransmitters.

Perhaps it could be worth studying whether the concept of interface, well defined in computer science, could be rethought in this context so as to facilitate the modular composition of stochastic models.

## Authors contributions

Each of the co-authors contributed equally to this project. More precisely, the first and the last author brought their competence in process calculi, and the second and the third their knowledge of neurophysiology, while they carried on jointly the design of the experiment, the modelization of the synapse and the analysis of the results.

## References

*Methods in Neuronal Modeling*. Koch C, Segev I edition. Cambridge: The MIT Press; 1998.*Computational Cell Biology*. Fall CP, Marland ES, Wagner JM, Tyson JJ edition. New York: Springer; 2002.**Rediscovering Biology - Molecular to Global Perspectives - Neurobiology**. [http://www.learner.org/channel/courses/biology/support/10_neuro.pdf]Frotscher M, Gundelfinger E, Jonas P, Neher E, Seeburg P:

**The most important recent advances in synsapse research from my point of view - and what remains to be done.***Cell Tissue Res*2006,**326:**203–204.Zucker RS, Kullmann DM, Schwartz TL:

**Release of Neurotransmitters.**In*From Molecules to Networks - An Introduction to Cellular and Molecular Neuroscience*. Edited by: Byrne JH, Roberts JL. San Diego: Elsevier - Academic Press; 2004:197–244.Sudhof TC:

**The synaptic vesicle cycle.***Annual Review of Neuroscience*2004,**27:**509–547.Schneggenburger R, Neher E:

**Intracellular calcium dependence of transmitter release rates at a fast central synapse.***Nature*2000,**46:**889–893.Schneggenburger R, Neher E:

**Presynaptic calcium and control of vesicle fusion.***Current Opinion in Neurobiology*2005,**15:**266–274.Meinrenken CJ, Borst JGG, Sakmann B:

**The Hodgkin-Huxley-Katz Prize Lecture. Local routes revisited the space and time dependence of the Ca signal for phasic transmitter release at the rat calyx of Held.**In*The Journal of Physiology*. in press;Smolen PD, Baxter DA, Byrne JH:

**Mathematical Modeling and Analysis of Intracellular Signaling Pathways.**In*From Molecules to Networks - An Introduction to Cellular and Molecular Neuroscience*. Edited by: Byrne JH, Roberts JL. Elsevier - Academic Press; 2004:391–429.Wilkinson DJ:

*Stochastic Modelling for System Biology*. London: Chapman and Hall - CRC Press; 2006.Gillespie DT:

**Exact stochastic simulation of coupled chemical reactions.***The Journal of Physical Chemistry*1977,**81:**2340–2361.Savtchenko LP, Rusakov DA:

**The optimal height of the synaptic cleft.***PNAS*2007,**104:**1823–1828.Graham B, Wong AYC, Forsythe ID:

**A computational model of synaptic transmission at the calyx of Held.***Neurocomputing*2001,**38–40:**37–42.Wong AY, Graham BP, Billups B, Forsythe ID:

**Distinguishing between presynaptic and postsynaptic mechanisms of short-term depression during action potentials trains.***J Neurosci*2003,**23:**4868–4877.Destexhe A, Mainen ZF, Sejnowski TJ:

**Synthesis of models for Excitable membrane, synaptic transmission and neuromodulation using a common kinetic formulation.***The Journal of Computational Neuro science*1994,**1:**195–231.Holmes WR:

**Calcium Signaling in Dendritic Spines.**In*Modeling in the Neuroscience - From Biological Systems to Neuromimetic Robotics*. 2nd edition edition. Edited by: Reeke GN, Poznanski RR, Lindsay KA, Rosenberg JR, Sporns O. NewYork: CRC Press; 2005:25–60.Hillston J:

*A Compositional Approach to Performance Modelling Cambridge:*. Cambridge University Press; 1996.Regev A, Shapiro E:

**Cellular Abstractions: Cells as Computation.***Nature*2002,**419:**343.Kitano H:

**Systems Biology: a brief overview.***Science*2002,**295**(5560):1662–1664.Cassman M, Arkin A, Doyle F, Katagiri F, Lauffenburger D, Stokes C:

*System Biology - International Research and Development*. The Netherlands: Springer; 2007.Nagasaki M, Onami S, Miyano S, Kitano H:

**Bio-calculus: its concept and molecular interaction.***Genome Informatics*1999,**10:**133–143.Priami C, Regev A, Shapiro E, Silvermann W:

**Application of a stochastic name-passing calculus to representation and simulation of molecular processes.***Theoretical Computer Science*2004,**325:**141–167.Regev A, Panina E, Silverman W, Cardelli L, Shapiro E:

**Bioambients: An abstraction for biological compartements.***Theoretical Computer Science*2004,**325:**141–167.Cardelli L:

**Brane calculi-interactions of biological membranes.**In*Proceedings of Computational Methods in Systems Biology*. Edited by: Vincent V, Schachter V. Paris: Springer; 2004:257–280. LNCS 3082Priami C:

**Stochastic π-calculus.***The Computer Journal*1995,**36**(6):578–589.Lecca P, Priami C, Quaglia P, Rossi B, Laudanna C, Costantin G:

**A stochastic process algebra approach to simulation of autoreactive lymphocyte recruitment.***SIMULATION: Trans. of the society for modelling and simulation international*2004,**80**(4):273–288.Calder M, Gilmore S, Hillston J:

**Modelling the Influence of RKIP on the ERK Signalling Pathway Using the Stochastic Process Algebra PEPA.***Transactions on Computational System Biology VII*2006, 1–23. LNBI 4230Milner R:

*Communicating and Mobile Systems: The π-Calculus*. Cambridge University Press; 1999.Phillips A, Cardelli L:

**Efficient, Correct Simulation of Biological Processes in the Stochastic Pi-calculus.**In*Proceedings of Computational Methods in Systems Biology*. Edited by: Calder M, Gilmore S. Edinburgh: Springer; 2007:184–199. LNCS 4695Kierzek AM:

**Stocks: Stochastic kinetic simulations of biochemical system with Gillespie algorithm.***Bioinformatics*2002,**18:**470–481.Schneggenburger R, Forsythe ID:

**The calxy of Held.***Cell Tissue Res*2006,**326:**311–337.Bollmann JH, Sakmann B:

**Control of synaptic strength and timing by the release-site***Ca*^{2+}**signal.***Nature Neuroscience*2005,**8:**426–434.Destexhe A, Mainen ZF, J ST:

**Kinetic Models of Synaptic Transmission.**In*Methods in Neuronal Modeling*. Edited by: Koch C, Segev I. Cambridge: The MIT Press; 1998:1–25.Clements JD:

**Transmitter timecourse in the synaptic cleft: its role in central synaptic function.***TINS*1996,**19:**163–171.Neher E:

**A comparison between exocytic control mechanisms in adrenal chromaffin cells and a glutamatergic synapse.***Pflugers Arch - Eur J Physiol*2006,**453:**261–268.

## Acknowledgements

This work has been partially supported by the MIUR project *Bisca* and by the Microsoft Research-Trento University Centre for Computational and Systems Biology.

This article has been published as part of *BMC Bioinformatics* Volume 9 Supplement 4, 2008: A Semantic Web for Bioinformatics: Goals, Tools, Systems, Applications. The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/9?issue=S4.

## Author information

### Authors and Affiliations

### Corresponding author

## Additional information

### Competing interests

The authors declare that they have no competing interests.

## Rights and permissions

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## About this article

### Cite this article

Bracciali, A., Brunelli, M., Cataldo, E. *et al.* Stochastic models for the in silico simulation of synaptic processes.
*BMC Bioinformatics* **9**
(Suppl 4), S7 (2008). https://doi.org/10.1186/1471-2105-9-S4-S7

Published:

DOI: https://doi.org/10.1186/1471-2105-9-S4-S7