Keywords

1 Introduction

The activity of neurons [16], including those taking place in the presynaptic bouton of the neuron, may be described theoretically with differential equations [1, 2, 4] or, more specifically, with partial differential equations (PDE) [3, 11, 13]. Numerical methods have been commonly used to solve PDE [9, 10, 12, 13], which may be a proof that scientists are constantly interested in this phenomenon and, on the other hand, that results are variable and depend strongly on the biological assumptions and on the way of modeling. To obtain the solution, the appropriate design of the mesh of the studied object (in this paper: of a presynaptic bouton of the neuron) is necessary, which may be a demanding task [19, 20].

Some experiments in the cited papers have already been performed in order to answer the question: how the neurotransmitter (NT) mediates the process of conducting nerve impulses, which is connected with the distribution of NT in the synapse and with its exocytosis. This paper is intended to expand that knowledge by using more complicated mathematical model in three dimensions together with a geometric model which is more complex than it was in some previous works [8, 13]. There are also examples of the description of a realistic geometric model [7], but without performing simulations. However, the NT flow in the realistically shaped model of the presynaptic bouton affected by the presence of a mitochondrion inside it has not been studied yet. Therefore, the objective of this paper was to examine the process of NT flow in a realistic model of the presynaptic bouton with a mitochondrion partly occupying its volume.

2 Theoretical Foundations

In this section the model based on partial differential equations is presented briefly. Such approach allows us to study the dynamics of the transport both in time and in spatial aspect. The model is nonlinear, diffusive-like - see the paper [3] for details. There are a few assumptions about the model. The bouton location is a bounded domain \(\varOmega \subset \mathbf {R}^3\). The boundary (the membrane) of the bouton is denoted as \(\partial \varOmega \). The total amount of neurotransmitter in the bouton is increasing when new vesicles are supplied inside the bouton and decreasing when their contents are released through the bouton membrane. The proper subset of the bouton, \(\varOmega _1\), is the domain of vesicles supply. The release site, \(\varGamma _d\), is a proper subset of the membrane. The function \(f:\varOmega \rightarrow \mathbf {R}\) models the synthesis of the vesicles that contain neurotransmitter. The value of this function is \(f(x) = \beta > 0\) on \(\varOmega _1\) and \(f(x) = 0\) on \(\varOmega \setminus \varOmega _1\). The neurotransmitter flow in the presynaptic bouton was modeled with the following partial differential equation:

$$\begin{aligned} \varrho _t(x,t) = a\varDelta \varrho (x,t) + f(x)(\bar{\varrho } - \varrho (x,t))^+ \quad \mathrm{in }\quad \varOmega \times (0,T), \end{aligned}$$
(1)

where \(\varrho (x,t)\) is the density of neurotransmitter at the point x and time t and a is the diffusion coefficient. The last term contains the function f(x), presented before, and the threshold NT density, above which synthesis dose not take place, denoted by \(\bar{\varrho }\). For any \(x\in \mathbf {R}\) the “\(x^+\)” symbol means max(0, x). The boundary conditions are

$$\begin{aligned} -\frac{\partial \varrho (x,t)}{\partial \nu } = 0 \quad \mathrm{on }\quad (\partial \varOmega \setminus \varGamma _d)\times (0,T), \end{aligned}$$
$$\begin{aligned} -\frac{\partial \varrho (x,t)}{\partial \nu } = \eta (t) \alpha \varrho (x,t)\quad \mathrm{on }\quad \varGamma _d\times (0,T), \end{aligned}$$
(2)

where \(\partial /\partial \nu \) is a directional derivative in the outer normal direction \(\nu \). The function \(\eta (t)\) depends on the time t and takes the value of 1 for \(t\in [t_n,t_n+\tau ]\) (with the action potential arriving at \(t_n\) and with the release time \(\tau \)), and \(\eta (t)=0\) otherwise [5, 8].

Multiplying (2) by a smooth function \(v:\varOmega \rightarrow \mathbf {R}\), and integrating by parts gives (see [3] for details):

$$\begin{aligned} \begin{aligned}&\int _{\varOmega }\varrho _t(x,t) v(x)\, dx + a\int _{\varOmega }\nabla \varrho (x,t)\cdot \nabla v(x)\, dx \\&\qquad +\alpha \eta (t)\int _{\varGamma _d}\varrho (x,t)v(x)\, dS =\beta \int _{\varOmega _1}(\bar{\varrho }-\varrho (x,t))^+v(x)\, dx, \end{aligned} \end{aligned}$$
(3)

provided that \(\varrho \) is sufficiently smooth.

Piecewise linear \(C^0\) tetrahedral finite elements are used in the numerical scheme. The unknown function \(\varrho \) is approximated by \(\varrho _h\), computed with the formula

$$\begin{aligned} \varrho _h (x,t) = \sum _{i=1}^K \varrho _i(t) v_i(x), \end{aligned}$$
(4)

in which the basis functions of the finite element space are denoted by \(v_i\). Using the substitution \(v=v_j\) in (2) yields

$$\begin{aligned} \begin{aligned}&\sum _{i=1}^K\varrho _i'(t)\int _{\varOmega }v_i(x) v_j(x)\, dx + a\sum _{i=1}^K\varrho _i(t)\int _{\varOmega }\nabla v_i(x)\cdot \nabla v_j(x)\, dx \\&\qquad + \alpha \eta (t)\sum _{i=1}^K\varrho _i(t)\int _{\varGamma _d}v_i(x)v_j(x)\, dS =\beta \int _{\varOmega _1}\left( \bar{\varrho }-\sum _{i=1}^K \varrho _i(t) v_i(x)\right) ^+v_j(x)\, dx. \end{aligned} \end{aligned}$$
(5)

When we denote \(\varrho _h(t) = [\varrho _1(t), \ldots , \varrho _K(t)]\), then the mass matrix \(\mathcal {G} = \{ g_{ij} \}_{i,j=1}^K\), the stiffness matrix \(\mathcal {A}_1 = \{ a^1_{ij} \}_{i,j=1}^K\), and the release matrix \(\mathcal {A}_2 = \{ a^2_{ij} \}_{i,j=1}^K\) are given by

$$\begin{aligned}&g_{ij} = \int _{\varOmega }v_i(x) v_j(x)\, dx, \\&a^1_{ij} = \int _{\varOmega }\nabla v_i(x)\cdot \nabla v_j(x)\, dx, \nonumber \\&a^2_{ij} = \int _{\varGamma _d}v_i(x)v_j(x)\, dS. \nonumber \end{aligned}$$
(6)

Let \(\mathcal {A}(t) = \mathcal {A}_1 + \eta (t)\mathcal {A}_2\). Let us also introduce the nonlinear operator \(\mathcal {B}:\mathbf {R}^K\rightarrow \mathbf {R}^K\) in the following way

$$\begin{aligned}&\mathcal {B}(\varrho _1,\ldots ,\varrho _K)_j = \beta \int _{\varOmega _1}\left( \bar{\varrho }-\sum _{i=1}^K \varrho _i v_i(x)\right) ^+v_j(x)\, dx. \end{aligned}$$
(7)

The Eq. (3) may be rewritten as

$$ \mathcal {G}\varrho _h'(t) + \mathcal {A}(t)\varrho _h(t) = \mathcal {B}(\varrho _h(t)). $$

The time step is defined as \(t_k = k \varDelta t\) where \(\varDelta t\) is at least several times smaller than \(\tau \). If we assume that \(\varrho _h^k\) is approximated by \(\varrho _h(t_k)\) we can use the Crank–Nicolson approximative scheme

$$\begin{aligned}&\mathcal {G}\frac{\varrho _h^{k}-\varrho _h^{k-1}}{\varDelta t} + \frac{1}{2}\mathcal {A}(t_k)\varrho _h^{k} +\frac{1}{2}\mathcal {A}(t_{k-1})\varrho _h^{k-1} = \frac{1}{2}\mathcal {B}(\varrho _h^{k}) + \frac{1}{2}\mathcal {B}(\varrho _h^{k-1}), \end{aligned}$$
(8)

valid for \(k = 1, 2, \ldots \). In the finite element basis \(\{ v_k \}_{k=1}^K\), \(\varrho _0\) is approximated by \(\varrho ^0_h\). The scheme must be solved for \(\varrho _h^k\). The problem of the nonlinearity in \(\mathcal {B}\) is dealt with by the iterative scheme:

$$\begin{aligned}&\varrho _h^{k(0)} = \varrho _h^{k-1},\\&\mathcal {G}\frac{\varrho _h^{k(m+1)}-\varrho _h^{k-1}}{\varDelta t} + \frac{1}{2}\mathcal {A}(t_k)\varrho _h^{k(m+1)} +\frac{1}{2}\mathcal {A}(t_{k-1})\varrho _h^{k-1} = \frac{1}{2}\mathcal {B}(\varrho _h^{k(m)}) + \frac{1}{2}\mathcal {B}(\varrho _h^{k-1}), \nonumber \end{aligned}$$
(9)

In each step of the iteration the linear system (9) is solved. The iterative procedure is stopped for the value m for which the correction of the solution due to one step of the scheme is sufficiently small.

Note that the non-linear term is calculated at the previous iteration (m) and the linear terms at iterations (m + 1). Such scheme reduces to solving the linear system in each step of the iteration. It is an alternative to the Newton method, and it has the advantage of not needing to cope with the non-differentiable nonlinearity introduced with the last term (taking the positive part) of the equation.

From physiological point of view, the time step for the Eq. (9) should be less than the time scale of the modeled phenomenon. We do not have the proof of the scheme stability but the numerical experiments show that the iterations of the fixed point scheme with respect to index m converge, and the limit is the solution of the Crank–Nicolson scheme. This scheme is known to be unconditionally stable. Still, to avoid numerical artifacts, the time-step needs to be sufficiently small. Our numerical experiments show that no numerical artifacts are seen in the obtained solution for \(\varDelta t = 0.0001\,{\upmu }\)s and for \(\varDelta t = 0.00005\,{\upmu }\)s.

For the numerical simulations, the bouton is modeled as a bounded polyhedral domain \(\varOmega \subset \mathbf {R}^3\). The boundary of the bouton, \(\partial \varOmega \), is represented by a set of flat polygons. The proper polyhedral subset of the bouton, \(\varOmega _1\), is the domain of vesicles supply. The release site \(\varGamma _d\) is modeled as a finite sum of flat polygons.

Various geometrical models of the presynaptic bouton have been used by the authors to simulate the process of conducting nerve impulses so far. The natural reference point for assessing their quality is to compare them to the images of real boutons which have been thoroughly examined, for example, several years ago [17, 18]. One of the models was based of the skeleton made up of two concentric spheres, the outer one having the diameter of 2–3 \({\upmu }\)m and the inner one free from synaptic vesicles, and also with numerous release sites [13]. The other structure consisting of two concentric globes and a single release site was filled with a tetrahedral mesh and it has been used to examine the distribution of neurotransmitter inside the bouton [6]. A very similar structure has been utilized to find the connection between the number and location of synthesis domains and the amount of neurotransmitter in different locations in the bouton [8]. Another, more complicated, surface model intended for the discrete simulation of vesicle movement consisted of a realistically shaped bouton and a mitochondrion partly occupying its space [7].

The amplitudes of evoked excitatory junctional potentials and the estimated number of the synaptic vesicles released during exocytosis were examined in [13] with the use of the model with empty central region. The study discussed, among others, the way the synaptic depression was reflected in the results of the simulations.

The model made up from two concentric spheres filled with a tetrahedral mesh was used to demonstrate synaptic depression [6]. The volume of the model was about 17,16 \({\upmu }\)m\(^3\) and the surface about 32,17 \({\upmu }\)m\(^2\). The diffusion coefficient was 3 \({\upmu }\)m\(^2\)/s, synthesis rate and the exocytosis rate around 21 \({\upmu }\)m\(^3\)/s. The results confirmed that during 0.5 s of 40 Hz stimulation, with the assumed values of parameters, the amount of neurotransmitter in the presynaptic bouton decreased, though only by the fraction of about 2\(\%\).

The next studies with the globe model [8] were performed to examine the influence of the location and number of synthesis zones on the amount of neurotransmitter in particular locations of the bouton. The simulation parameters were similar as in [6] but the number of synthesis zones was 1 or 2 (however, with the same total volume) and the location was also variable. The chief conclusion was that the closer the synthesis (or supply) domain to the release site, the faster the bouton becomes depleted of neurotransmitter.

3 The Realistic Model of the Bouton with Mitochondrion

The geometric aspects of the models of the bouton used so far, mentioned in the previous section, were extremely simplified. In this paper, the simulations have been conducted on the basis of the realistic model of the bouton in which the mitochondrion is modeled as well.

3.1 The Input Surface Mesh

The structure used to model the presynaptic bouton in this paper was based on the real shape of such a bouton and it was composed of geometric shapes modified so that it resembled the original. Therefore, the outer of the bouton was modeled as a moderately modified sphere, whereas a part of the mitochondrion inside it was a highly distorted cuboid, as it has been presented in [7]. The parameters of the surface mesh are described therein in detail. It should be stressed that the mesh described in this paper was designed to shorten the computing time of simulations of neurotransmitter flow, and therefore contains, among others, relatively large elements. The number of tetrahedra in the mesh was less then \(10^5\) and the number of the surface triangles (faces) was less than \(5\times 10^4\).

3.2 Three-Dimensional Model

The parameters of the input surface mesh were as follows. The total surface of the bouton was \(S\,{\approx }\,6.7811\,\) \({\upmu }\)m\(^2\), and the area of the active zone (of the release site) was \(S_{AZ}\,{\approx }\,0.2402\,{\upmu }\)m\(^2\). The surface mesh is presented in Fig. 1.

Fig. 1.
figure 1

The surface mesh of the presynaptic bouton - the view with the top cut off. The bouton would fit in a 1.33 \(\times \) 1.16 \(\times \) 1.23 \({\upmu }\)m cuboid. Skeleton - left, partly colored - middle, colored surface - right. The meaning of the colors is: red(-4) - mitochondrion, green(-3) - artificial cut, blue(-2) - release site, yellow(-1) - the rest. (Color figure online)

Full size image

The tetrahedral mesh was generated with the TetGen program [14, 15]. The result is depicted in Fig. 2.

The input parameters for TetGen were chosen in order to minimize the number of mesh elements. As a result of this, the bouton three-dimensional mesh contained 77418 tetrahedra with the total volume \(V\,{\approx }\,0.9029\,{\upmu }\)m\(^3\). The volume of the part of the bouton located in the direct neighborhood of the mitochondrion (where the endings of the microtubules are sometimes found), assumed as the theoretical NT supply zone (“synthesis domain”) was \(V_s\,{\approx }\,0.0198\,{\upmu }\)m\(^3\). The quality of the mesh was assessed by computing several measures for each tetrahedra and by taking, for each measure, its maximal value i.e. the value for the worst mesh element. The values are collected in Table 1. The relatively high values of the mesh quality parameters may suggest a cautious approach to the analysis of the results. However, further stages of the experiment revealed that the mesh proved sufficient for performing the planned simulation.

Table 1. The parameters of the quality of three dimensional mesh of the presynaptic bouton with a mitochondrion.
Full size table
Fig. 2.
figure 2

The cross section of the tetrahedral mesh of the bouton used in calculations. The bouton would fit in a 1.33 \(\times \) 1.16 \(\times \) 1.23 \({\upmu }\)m cuboid. The meaning of the colors: red(-4) - inner surface of the mitochondrion, green(-3) - artificial cut, blue(-2) - release site, yellow(-1) - the rest of the bouton surface, purple - the elements of the tetrahedral mesh. (Color figure online)

Full size image

The initial density of the neurotransmitter was calculated by using the formula

$$\begin{aligned} {} {\varrho }(t=0)={\varrho }_0=Ae^{-br^2} \text {[vesicles/}{\upmu }\text {m}^3] \end{aligned}$$
(10)

where \(A=300\)[vesicles/\({\upmu }\)m\(^3\)] is the theoretical maximal value of the function, \(b=0.28[1/{\upmu }\)m\(^2\)], and \(r[{\upmu }\)m] is the distance from the “center” of the bouton i.e. from the center of the minimal cuboid containing the model.

4 Simulation Results

The simulation time was 0.1 s and during that time there were 4 impulses. The results are depicted in the following figures: Fig. 3 and Fig. 4, and in Table 2.

Table 2. The amount of neurotransmitter in vesicles during 4 consecutive nerve impulses. The columns contain: Time - simulation time, Total - the total amount of NT in the bouton, Rel - the amount of NT released in a single time step. The NT amount unit is the amount of NT in a single synaptic vesicle.
Full size table

The program was designed primarily for validating the accuracy of the calculations; therefore the algorithm was implemented as a single threaded code in Python, with numeric and graphic modules included, and the calculations were relatively slow. The approximate speed was 30–40 iterations per hour. Therefore one simulation run lasted about 2–3 days.

The total amount of neurotransmitter in the presynaptic bouton calculated during simulation was almost constant. The relative decrease of its value did not exceed 1\(\%\). Those values refer to the situation when the neuron is not very intensively stimulated, and synaptic depression is not likely to occur. From the analysis of Fig. 3 it may be concluded that the activity of the release site is very moderate. However, closer examination of Fig. 4 reveals that the spatial distribution of neurotransmitter does change during stimulation. In the region directly adjacent to the release site the amount of neurotransmitter drops a little, and also in its vicinity a slight decrease may be noticed, which is visible in the last picture in Fig. 4. The amount of neurotransmitter increased between stimuli, though at a very low speed, which may be attributed to the low value of exocytosis rate. The changes in these parameters can result in neurotransmitter amount sinking rapidly during exocytosis, thus leading to synaptic depression.

To verify the reliability of the results, two control simulations were run. The first one, with two times larger mesh, confirmed that the relative differences in the results did not exceed 0.14\(\%\). For the second one, with two times larger mesh (as before) and with halved time step, the results were almost the same; the relative difference between two control runs did not exceed \(10^{-10}\).

Fig. 3.
figure 3

The spatial distribution of neurotransmitter density ND(x, y, z) in the presynaptic bouton after the arrival of 3 consecutive action potentials. The total time was 0.1 s and the impulses, presented in the rows, arrived at time intervals 0.0123 s–0.0127 s, 0.0373 s–0.0377 s, and 0.0623 s–0.0627 s. In the last row the first (0.0 s) and the last (0.1 s) moment of simulation is depicted. The tick marks in all axes are placed from −0.6 \({\upmu }\)m to +0.6 \({\upmu }\)m, at 0.2 \({\upmu }\)m intervals. The meaning of colors, corresponding to the range of ND (0\(\div \)350 vesicles/\({\upmu }\)m\(^3\)) is explained by the color bar at the top. (Color figure online)

Full size image
Fig. 4.
figure 4

The radial distribution (\(\varrho (r)\), where \(r=\sqrt{x^2+y^2+z^2}\)) of NT density in the presynaptic bouton after 3 consecutive action potentials. Three impulses are in the first 3 rows, and in the last row the plots for the beginning (left) and for the end (right) of the simulation are drawn. The total time was 0.1 s and the impulses time intervals 0.0123 s–0.0127 s, 0.0373 s–0.0377 s, and 0.0623 s–0.0627 s. The tick marks in the r axis are in the range 0 \(\div \)0.75 \({\upmu }\)m, placed each 0.1 \({\upmu }\)m, and in the \(\varrho \) axis they are placed from 0 to 350 vesicles/\({\upmu }\)m\(^3\) at the intervals of 100 vesicles/\({\upmu }\)m\(^3\). The meaning of colors, reflecting the z coordinate (−0.75 \(\div \)+0.75 \({\upmu }\)m) is explained by the color bar at the top. (Color figure online)

Full size image

The process described throughout this paper is similar to which has been found before [1] where, in Fig. 4, one can notice the almost constant number of released vesicles in every time step at the end of the simulation. Taking into account the fact that the authors of the cited paper assumed that the total number of synaptic vesicles in the bouton was more than \(10^4\), and that the number of synaptic vesicles released in each time step was about 40, we may conclude that if with \(M\,\approx \,250\) vesicles in the bouton (the value assumed in this paper) the amount of neurotransmitter released during one time step is approximately 0.5, our results are similar to those found in literature; the proportions of the released amount to the total vesicle pool are of the same order of magnitude. In another study [2] the equilibrium state has been achieved, though at a lower level, indicating synaptic depression.

5 Concluding Remarks

The model of realistically shaped presynaptic bouton with a mitochondrion partly blocking its volume, presented in this paper, proved its ability to be used in simulation of synaptic depression. It should be stressed that the values of the parameters chosen for the initial tests of the proposed structure refer to the situation of a very regular activity, not threatened by depression. However, in silico tests may reveal the changes in distribution of neurotransmitter in a presynaptic bouton during a very frequent stimulation, thus allowing us to study depression in detail, the more so because the results of the experiments found in literature, whether or not a strong synaptic depression was detected, confirm that our simulation results reflect the real processes taking place in the presynaptic bouton of a stimulated neuron.