Relations between Geometric Parameters and Numerical Simulation Accuracy in Modeling Signal Transmission in the Presynaptic Bouton

Featured Application

The methodology described in this paper may be used, and will be used by the authors of this and related works, to (1) examine the influence of various factors, including the shape of the bouton and its neurotransmitter synthesis domains as well as the release sites and various organelles, on the course of the nerve impulse. And, vice versa, (2) to detect some abnormalities in the composition o the synapse basing on the shape of its electrical output.

Abstract

In order to examine nerve impulses by means of simulation methodology, the models of all parts of nervous system, well suited for numerical modeling, are needed. In this paper the problem of setting up such a model, namely, that of a presynaptic bouton, is addressed. Simulation of the neurotransmitter flow inside the presynaptic bouton is performed. The transport is modeled with a partial differential equation with an additional nonlinear term. Two ways of modeling the bouton are applied. One of them let reflect a complex shape of the bouton and of some inner organelles. The influence of the generated mesh quality on the accuracy of numerical simulations is studied by comparing the released amount of neurotransmitter. The only mesh that produced diminished output was the worst one. The conclusion is that even slightly optimized tetrahedral mesh is suitable for calculations.

1. Introduction

Dynamical processes in the neuron [1] include, among others, those concerning the release of neurotransmitter (NT) into the synaptic cleft. The release occurs via calcium channels [2,3] and the amount of the released NT is measured indirectly by the postsynaptic potential. In Figure 1A in the cited paper [3] the connection between NT amount and postsynaptic current is illustrated; in Figure 3C in other work [4] the shape of the impulse in genetically different animals is presented; diminished amplitude means impaired genotype. The above mentioned dynamical processes are often described by using differential equations [5,6,7], including—for a long time—partial ones [8,9,10]. The equation of this type can not, usually, be solved analytically, and therefore numerical methods have to be used [10,11,12,13]. In order to perform the calculations quickly and accurately, the appropriately designed mesh, in the nodes of which the calculations will be made, should be created. In biological sciences the boundary conditions are generated by the modeled organelles in biological structure of the cell. For example, what the model should take into account is the work of the active zone of the presynaptic bouton during exocytosis. It is one of the crucial mechanisms in transmitting nerve impulses and it is accomplished by means of $C{a}_{V}2$ channel (one of three kinds of similar channels: $C{a}_{V}1$, $C{a}_{V}2$ and $C{a}_{V}3$) [4]. Therefore, the design of a good mesh is a difficult problem connected with the mesh geometrical properties [14,15]. Tetrahedral meshes are also studied in the context of automatic adjustment [16]. Finite Elements Method, on the other hand, requires a high quality mesh [17]. This means that a single cell should be as regular as possible—an equilateral triangle in two dimensions, and in the 3D case—a regular tetrahedron. Generally, it is impossible to reach this goal in practice because of the geometric properties of the modeled object. Additionally, regular tetrahedra can not be used as elements of a spatial mesh. In reality, such tetrahedra differ in shape and size; also, the regions joining different parts of the model are usually covered with elements of worse quality.

In this paper the problem of dependence of numerics precision on the quality of the spatial data structure (a tetrahedral mesh) for numerical calculations for processes in the neuron axon terminal is studied. Various geometrical presynaptic bouton models are discussed. They are compared by using several different parameters of the quality of the mesh, one of them formerly introduced [18]. The differences in the results of simulations for different meshes, with other parameters kept unchanged, are analyzed. The adequate measure of the differences is the amount of neurotransmitter released to the synaptic cleft.

This paper is a continuation of the studies described in [6,8,11,12,18,19,20,21,22].

2. Model

Since the transport of neurotransmitter (NT) in the presynaptic bouton of the neuron has to be considered both in spatial and in temporal aspect, the partial differential equation model has been introduced, derived from the diffusion equation with nonlinear modifications. The more detailed description of the model is presented in [8]. Main aspects of the model are given below.

2.1. Diffusive PDE Model

When setting up the model, the following was assumed:

The presynaptic bouton occupies a bounded subset $\Omega$ of three-dimensional space. The NT supply takes places in the bouton subdomain ${\Omega }_1$ (supply zone), thereby increasing the total amount of NT in the bouton. On the other hand, the release through the subset ${\Gamma }_d$ of the bouton membrane $\partial \Omega$ decreases that amount. The NT flow in the presynaptic bouton, basing on the similar study [22], was modeled with the following partial differential Equation (1):

$$\begin{array}{ccc}\hfill & & {\varrho }_t(x,t)=a\Delta \varrho (x,t)+f\left(x\right){(\overline{\varrho }-\varrho (x,t))}^{+}\mathrm{in}\Omega \times (0,T),\hfill \end{array}$$

(1)

where $\varrho (x,t)$ is the density of NT at time t and at the spatial location x and a is the diffusion coefficient. The appropriate initial conditions had to be set and, similarly, the boundary conditions [22] should be modeled by the Equation (2)

$$\begin{array}{ccc}\hfill & & -\frac{\partial \varrho (x,t)}{\partial \nu }=\eta \left(t\right)\alpha \varrho (x,t)\mathrm{on}{\Gamma }_d\times (0,T)\hfill \end{array}$$

(2)

where $\alpha$ is the permeability coefficient and $\partial /\partial \nu$ is a directional derivative in the outer normal direction $\nu$. The function $\eta \left(t\right)$ (depending on time t) 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 \left(t\right)=0$ otherwise [19]. In such a way, boundary conditions are a simple model of exocytosis through a $C{a}_{V}2$ type voltage-gated channel.

When the Equation (1) is multiplied by a smooth function $v:\Omega \to \mathbb{R}$, and after integrating by parts [22], we obtain

$$\begin{array}{cc}\hfill & {\int }_{\Omega }{\varrho }_t(x,t)v\left(x\right)dx+a{\int }_{\Omega }\nabla \varrho (x,t)·\nabla v\left(x\right)dx\hfill \\ \hfill & +\alpha \eta \left(t\right){\int }_{{\Gamma }_d}\varrho (x,t)v\left(x\right)dS=\beta {\int }_{{\Omega }_1}{(\overline{\varrho }-\varrho (x,t))}^{+}v\left(x\right)dx.\hfill \end{array}$$

(3)

which is valid when $\varrho$ is sufficiently smooth.

2.2. 3D Simulations—The Finite Elements Method

Piecewise linear $C^0$ finite elements (tetrahedra) are used to perform the calculations according to the numerical scheme [22]. The function ${\varrho }_h$, which approximates $\varrho$, is calculated as

$${\varrho }_h(x,t)=\sum _{i=1}^K{\varrho }_i\left(t\right)v_i\left(x\right),$$

where $v_i$ are the basis functions of the finite elements space. When the following matrices are introduced: The mass matrix $\mathcal{G}={\left\{g_{ij}\right\}}_{i,j=1}^K$, the stiffness matrix ${\mathcal{A}}_1={\left\{a_{ij}^1\right\}}_{i,j=1}^K$, and the release matrix ${\mathcal{A}}_2={\left\{a_{ij}^2\right\}}_{i,j=1}^K$, and when the additional assumption is made that $\mathcal{A}\left(t\right)={\mathcal{A}}_1+\eta \left(t\right){\mathcal{A}}_2$, and the nonlinear operator $\mathcal{B}:{\mathbb{R}}^K\to {\mathbb{R}}^K$ is applied according to [22]:

$$\begin{array}{ccc}& & \mathcal{B}{({\varrho }_1,\dots ,{\varrho }_K)}_j=\beta {\int }_{{\Omega }_1}{\left(\overline{\varrho }-\sum _{i=1}^K{\varrho }_i{v}_i\left(x\right)\right)}^{+}{v}_j\left(x\right)dx.\hfill \end{array}$$

(4)

the numerical scheme may be rewritten as

$$\mathcal{G}{\varrho }_h^{\prime }\left(t\right)+\mathcal{A}\left(t\right){\varrho }_h\left(t\right)=\mathcal{B}\left({\varrho }_h\left(t\right)\right).$$

The time step is defined as $t_k=k\Delta t$ where $\Delta t$ is at least several times smaller than $\tau$. After the approximation of ${\varrho }_h^k$ by ${\varrho }_h\left(t_k\right)$, the Crank–Nicholson approximative scheme can be used

$$\begin{array}{ccc}& & \mathcal{G}\frac{{\varrho }_h^k-{\varrho }_h^{k-1}}{\Delta t}+\frac{1}{2}\mathcal{A}\left(t_k\right){\varrho }_h^k+\frac{1}{2}\mathcal{A}\left(t_{k-1}\right){\varrho }_h^{k-1}=\frac{1}{2}\mathcal{B}\left({\varrho }_h^k\right)+\frac{1}{2}\mathcal{B}\left({\varrho }_h^{k-1}\right),\hfill \end{array}$$

(5)

valid for $k=1,2,\dots$. In the finite element basis ${\left\{v_k\right\}}_{k=1}^K$, ${\varrho }_0$ is approximated by ${\varrho }_h^0$. 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{array}{c}{\varrho }_h^{k\left(0\right)}={\varrho }_h^{k-1},\hfill \\ \mathcal{G}\frac{{\varrho }_h^{k(m+1)}-{\varrho }_h^{k-1}}{\Delta t}+\frac{1}{2}\mathcal{A}\left(t_k\right){\varrho }_h^{k(m+1)}+\frac{1}{2}\mathcal{A}\left(t_{k-1}\right){\varrho }_h^{k-1}=\frac{1}{2}\mathcal{B}\left({\varrho }_h^{k\left(m\right)}\right)+\frac{1}{2}\mathcal{B}\left({\varrho }_h^{k-1}\right),\hfill \end{array}$$

(6)

in which the inner iteration step is denoted by the index m. In such a way, the scheme gives a linear system of equations for each iteration. Iterating stops when $|{\varrho }_h^{k(m+1)}-{\varrho }_h^{k\left(m\right)}|$ is sufficiently small.

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

3. Generation of the Mesh

As it has already been mentioned [18], the problem of generating a mesh of a complex object is difficult. Firstly, a reasonable compromise must be found between a mesh with a very large number of small tetrahedra (increased computing time) and another one with a small number of relatively large elements (errors in simulating flow near tiny substructures). Secondly, the mesh elements should be as regular as possible [23].

3.1. Mesh Element Quality

The mesh quality parameters, considered below and calculated for each element, are as follows.

Assuming that S denotes the total surface, and $V$—the total volume of the tetrahedron, the following parameters have been used [18]: The first of the defined parameters was

$$SV=\frac{\sqrt{S}}{\sqrt[3]{V}}$$

The second parameter was the ratio of the longest edge and the radius of the inscribed sphere [24]:

$$ER=\frac{l_{max}}{r_{ins}},$$

ranging from $E{R}_{min}=2\sqrt{6}\approx 4.89898$, for a regular tetrahedron, to infinity.

The third one is the value applied in the TetGen program [25,26] i.e., the longest edge divided by the smallest height (EH)

$$EH=\frac{l_{max}}{h_{min}},$$

with the range similar as in the previous example i.e., from $E{H}_{min}=\sqrt{\frac{3}{2}}\approx 1.22474$ to infinity, following the same rule as before: “the larger the value, the worse the tetrahedron”.

The last two parameters are $MAXS$ and $MINS$, and they are computed as the ratios of the largest ($S_{max}$) or the smallest ($S_{min}$) face measure, respectively, to the total surface (S) of the tetrahedron.

The formula for the first of those two parameters is

$$MAXS=\frac{S_{max}}{S}$$

and its value ranges from 0.25 to 0.5 i.e.,

$$MAXS\in [0.25;0.5)$$

The formula for the second parameter is

$$MINS=\frac{S_{min}}{S}$$

and its value ranges from 0 to 0.25 i.e.,

$$MINS\in (0;0.25]$$

The parameters mentioned above (SV, ER, EH, MAXS and MINS) is positive and dimensionless. All of them, except MINS, have the smallest value for the regular tetrahedron. The MINS value behaves differently: the closer to zero, the worse the tetrahedron; the maximum value 0.25 is reached for a regular tetrahedron. Moreover, the value 0.25 of both the MAXS and the MINS does not necessarily mean that the tetrahedron is regular—it may be a sliver-type (see Figure 8 in [25]). Other irregular tetrahedron with the vertices at (−0.001,0,0), (0.001,0,0), (0,−0.001,100) and (0,0.001,100) is a wedge-type one (see Figure 8 in [25]).

The numerical approximation of a gradient of a given function depends strongly on the largest angle (planar angle for a triangle and spherical or dihedral angle for a tetrahedron) [27]. Thus, the condition number of the stiffness matrix should be kept to a minimum.

Theoretically, to solve a partial differential equation numerically, the above mentioned quality measures $Q_1$–$Q_5$ should be as close to the ideal ones as possible for each of the mesh elements. However, not all deviations from the ideal have strong negative influence on the calculations with FEM. Moreover, it turns out that tetrahedrons with faces far from optimal are still of a quality enough to be used for numerical analysis of the model [23].

3.2. Mesh Boundary

If the structure to model is complex in a sense that it contains many substructures with their own surfaces which consist of triangles of very different size and quality then the generated mesh must take into account these substructures and the joints between them, which constitutes additional problems during the mesh generation.

The 3D mesh should be of a good quality not only in its 3D aspect but also in its 2D aspect as well. This concerns the quality of the surface of the mesh and, especially, the parts belonging to distinct areas.

The quality of a 2D (triangular) mesh of the surface of the bouton which topologically remains similar to the 3D ball is assessed analogically to the quality of the 3D net.

The quality measures that correspond to the ones defined for the spatial meshes are as follows: the square perimeter divided by the square root of the surface (PS2), the longest edge divided by the radius of the inscribed sphere (in that case the 2D sphere is a circle) (ER2), the longest edge divided by the smallest height (EH2) or the function of the ratio of the largest and smallest edge length, respectively, to the total perimeter of the triangle (MAXE and MINE).

Similar as in the 3D case, the best (the smallest) attainable values of the above measures (except MINE) are:

$PS{2}_{min}=2x{27}^{1/4}\approx 4.55901$

$ER{2}_{min}=2x{3}^{1/2}\approx 3.4641$

$EH{2}_{min}=2/\left[3^{1/2}\right]\approx 1.1547$

$MAX{E}_{min}=1/3\approx 0.3333$

and the largest MINE is

$MIN{E}_{max}=1/3\approx 0.3333$

For the first three measures the upper bound is infinity whereas for the fourth one (MAXE) it is equal to 0.5, and the lower bound for the last parameter (MINE) is 0.

3.3. Mesh Inconsistency

One of the possible faults of the mesh results from the fact that its boundary is not a closed surface and the meshing program is trying to create elements (tetrahedra) outside the bouton (leaking) which may happen when the intended common vertices of the adjacent tetrahedra separate, creating the gap between the walls. Another example of “leaking”, termed as double edge, means that a cleft is created between the edges and the size of the cleft is less than one thousandth of the edge length. Intersection appears when the edge of one tetrahedron is intersecting with the face of another one; Some subset of the 3D mesh volume is then defined twice.

3.4. Mesh Quality

The measures of the mesh quality can by either the extreme values of the quality parameters or other parameters of the distribution of a parameter e.g., the maximum value of SV, ER, EH, or MINS, or minimum value of MAXS. These quantities, that are sometimes used to describe the mesh, and also the distributions themselves, are presented as tables and histograms.

4. Geometrical Models of the Bouton

The geometry of the bouton was modeled in three various ways, using three different objects.

Two of them—the geosphere and the ball (globe)—were extremely simplified models. The basic difference between geosphere surface mesh and globe surface mesh is that the geosphere surface consists only of triangles, nearly equally-shaped and nearly equilateral, whereas the globe surface contains two types of polygons: (1) equal isosceles (but not necessarily equilateral) triangles; and (2) trapezoids. The two above models are known to differ in some aspects. Two of them are most commonly mentioned. The first one is that with the globe mesh it is easier to create shapes with edges along meridians and circles of latitude. The second one is that with the geosphere mesh it is easier to make a relatively smooth surface with a smaller number of mesh elements.

The third model, so-called real model, reflected the real geometry of the bouton surface and accounted for the presence of the inner organelles—see [20] for details.

4.1. Geosphere—The Simplest Model of the Sphere

The two-dimensional geosphere surface mesh, often used, for example in [18], to depict the surface of a sphere, consists of triangles, which are almost equilateral. Their characteristics exceed the ideal values by 1.3% to 25.5% which constitutes a good reference point for comparison of the quality of 3D meshes. The analyzed surface (Figure 1, left) had 15,680 triangular faces and 7482 vertices.

The quality of the model of the geosphere surface is presented in

Table 1. Geosphere model surface mesh quality.

Name	Definition	Ideal	Min	Max	Max/Ideal
$PS2$	$P/\sqrt{S}$	$2\sqrt[4]{27}\approx 4.55901$	$4.55903$	$4.6104$	$1.011$
$ER2$	$E_{max}/R$	$2\sqrt{3}\approx 3.4641$	$3.4684$	$3.9245$	$1.133$
$EH2$	$E_{max}/H_{min}$	$2/\sqrt{3}\approx 1.1547$	$1.1575$	$1.4527$	$1.258$
$MAXE$	$E_{max}/P$	$1/3\approx 0.3333$	$0.3337$	$0.3702$	$1.111$
$MINE$	$E_{min}/P$	$1/3\approx 0.3333$	$0.2913$	$0.3325$	1.144 ${}^1$

¹ This value is calculated differently—as ideal/min.

.

The angles of the triangular mesh considered in Table 1 vary between the values of $48.53°$ and $71.98°$. The range of triangle area was from $S_{min}=0.001903{\mathrm{\mu m}}^2$ to $S_{max}=0.002051{\mathrm{\mu m}}^2$. Therefore, the ratio $S_{max}/S_{min}$ was $107.8\%$. The surface mesh was then submitted as the input data to the TetGen program which produced the tetrahedral mesh of the sphere (Figure 1, right).

The parameters of the sphere were as in

Table 2. Mesh quality of the geosphere 3D model.

Name	Definition	Ideal	Min	Max	Max/Ideal
$SV$	$\sqrt{S}/\sqrt[3]{V}$	$\sqrt[4]{3}\sqrt[6]{72}\approx 2.6843$	$2.6844$	$4.7147$	$1.756$
$ER$	$E_{max}/R$	$2\sqrt{6}\approx 4.8989$	$4.9301$	$37.4427$	$7.643$
$EH$	$E_{max}/H_{min}$	$\sqrt{3}/\sqrt{2}\approx 1.2247$	$1.2398$	$14.628$	$11.944$
$MAXS$	$S_{max}/S$	$1/4=0.2500$	$0.2503$	$0.4434$	$1.774$
$MINS$	$S_{min}/S$	$1/4=0.2500$	$0.0639$	$0.2498$	3.912 ${}^1$

¹ This value is calculated differently—as ideal/min.

.

These parameters are the reference point for the next calculations.

4.2. The Ball (Globe)—The Simplest Geometric Model of the Presynaptic Bouton

The globe is a sphere divided into spherical polygons by meridians and circles of latitude. A typical globe is a depiction of the Earth’s surface. For example, one may draw a globe with 24 meridians and 13 circles of altitude—or rather 11, if we exclude poles as degenerate circles. The simplest object used for modeling a presynaptic bouton [19] consisted of two such globes placed concentrically (see Figure 2).

The input parameters for the meshing software were as follows: maximum volume of a single mesh element (tetrahedron) ($V{T}_{max}$). Another two important parameters are the maximum aspect ratio ($E{H}_{max}$) which is not standardized but computed directly as the ratio of the longest edge and the shortest height of a tetrahedron. The last restriction is that of the minimum dihedral angle ($D{A}_{min}$). It should be stressed that the last restriction is the most strictly observed in the program.

The 2-dimensional mesh is generated in such a way that each spherical polygon is replaced by a flat polygon with vertices at the same points. The quality of the globe mesh was tested in the following subspace of the above parameters:

$$E{H}_{max}\in [4.0;10.0]$$

$$V{T}_{max}\in [0.0006;4.0]$$

$$D{A}_{min}\in [10°;19°]$$

The example results of these tests are presented in Figure 3:

• for the fixed values of $E{H}_{max}=5$ and $V{T}_{max}=0.001$ and for $D{A}_{min}$ ranging from $10°$ to $19°$;
• for the fixed values of $E{H}_{max}=5$ and $D{A}_{min}=13°$ and for $V{T}_{max}$ ranging from 0.0006 to 0.0015.

The best quality was achieved when the imposed restrictions were: $E{H}_{max}=10.0$, $D{A}_{min}=13°$ and $V{T}_{max}=0.0012$.

The input 2D mesh and the results of its tetrahedralization are presented in

Table 3. Mesh quality of the globe surface model.

Name	Definition	Ideal	Min	Max	Max/Ideal
$PS2$	$P/\sqrt{S}$	$2\sqrt[4]{27}\approx 4.55901$	$4.55910$	$6.3549$	$1.394$
$ER2$	$E_{max}/R$	$2\sqrt{3}\approx 3.4641$	$3.4792$	$9.6114$	$2.775$
$EH2$	$E_{max}/H_{min}$	$2/\sqrt{3}\approx 1.1547$	$1.1647$	$4.5749$	$3.962$
$MAXE$	$E_{max}/P$	$1/3\approx 0.3333$	$0.3348$	$0.4575$	$1.429$
$MINE$	$E_{min}/P$	$1/3\approx 0.3333$	$0.1617$	$0.3325$	2.061 ${}^1$

¹ This value is calculated differently—as ideal/min.

and

Table 4. Mesh quality of the globe 3D model.

Name	Definition	Ideal	Min	Max	Max/Ideal
$SV$	$\sqrt{S}/\sqrt[3]{V}$	$\sqrt[4]{3}\sqrt[6]{72}\approx 2.68433$	$2.68503$	$4.23589$	$1.578$
$ER$	$E_{max}/R$	$2\sqrt{6}\approx 4.89898$	$4.97976$	$28.12424$	$5.741$
$EH$	$E_{max}/H_{min}$	$\sqrt{3}/\sqrt{2}\approx 1.22474$	$1.26564$	$11.40720$	$9.314$
$MAXS$	$S_{max}/S$	$1/4=0.25$	$0.25037$	$0.44492$	$1.780$
$MINS$	$S_{min}/S$	$1/4=0.25$	$0.06135$	$0.24977$	4.075 ${}^1$

¹ This value is calculated differently—as ideal/min.

.

4.3. The Real Model with Inner Organelles

In the real model, the shape of the bouton surface is created on the basis of the microscope image [28]—see [20] for details concerning the way the model is generated. Furthermore, the inner organelles—mitochondrion and microtubules—are also modeled as the inner structures of the bouton—Figure 4. In the model the following areas have to be taken into consideration: the area of the synapse membrane in general, the zone of releasing the neurotransmitter particles, the surface of the mitochondrion which in our model is connected to the artificially (i.e., numerically) cut bouton, the surface of the microtubules and of microtubule terminals. The above surfaces are marked with different colors in the graphs—Figure 5.

The input mesh was made up of 41,807 vertices and 83610 triangles—Figure 6. Generation of the three-dimensional mesh in such a way that the triangular faces of the mesh tetrahedra constitute the two-dimensional mesh covering the surfaces of the inner structures, especially when small details are concerned, was one of the crucial problems (Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, in particular Figure 5 and Figure 9, see also [17]). The characteristics of the input mesh are presented in

Table 5. Surface mesh quality of the real 3D model.

Name	Definition	Ideal	Min	Max	Max/Ideal
$PS2$	$P/\sqrt{S}$	$2\sqrt[4]{27}\approx 4.55901$	$4.55902$	$21.12052$	$4.633$
$ER2$	$E_{max}/R$	$2\sqrt{3}\approx 3.4641$	$3.4684$	$111.47228$	$32.179$
$EH2$	$E_{max}/H_{min}$	$2/\sqrt{3}\approx 1.1547$	$1.15740$	$55.71270$	$48.249$
$MAXE$	$E_{max}/P$	$1/3\approx 0.33333$	$0.33372$	$0.49979$	$1.499$
$MINE$	$E_{min}/P$	$1/3\approx 0.3333$	$0.04002$	$0.33296$	8.329 ${}^1$

¹ This value is calculated differently—as ideal/min.

.

The above input mesh was passed as input data to the TetGen program. The program generated output mesh with no optimization, thereby the result was deliberately of a very poor quality. The parameters of the output tetrahedral mesh are presented in

Table 6. Mesh quality of the real 3D model—the worst mesh.

Name	Definition	Ideal	Min	Max	Max/Ideal
$SV$	$\sqrt{S}/\sqrt[3]{V}$	$\sqrt[4]{3}\sqrt[6]{72}\approx 2.684$	$2.703$	$52.387$	$19.516$
$ER$	$E_{max}/R$	$2\sqrt{6}\approx 4.8990$	$5.7381$	$\mathrm{52,443.6}$	$\mathrm{107,049.6}$
$EH$	$E_{max}/H_{min}$	$\sqrt{3}/\sqrt{2}\approx 1.2247$	$1.5641$	$\mathrm{262,212.2}$	$\mathrm{214,095.4}$
$MAXS$	$S_{max}/S$	$1/4=0.25000$	$0.25002$	$0.499995$	$1.99998$
$MINS$	$S_{min}/S$	$1/4=0.25000$	$0.000009$	$0.24998$	28,444.1 ${}^1$

¹ This value is calculated differently—as ideal/min.

. The mesh had 57,060 nodes, 114,116 boundary faces and 194,064 tetrahedra.

The next step was, as it has been done previously, to generate the optimized mesh. Its quality was checked for different values of input TetGen parameters. Those parameters after preliminary tests were set to

$E{H}_{max}\in [70;100]$

$D{A}_{min}\in [1°;20°]$

$V{T}_{max}$ = 1000 $\mathsf{\mu }{\mathrm{m}}^3$

The best quality was achieved for $E{H}_{max}=100;D{A}_{min}=2°;V{T}_{max}=1000$.

The size of the mesh was 391,897 nodes, 4,141,116 faces (296,528 on the bouton surface) and 1,996,426 tetrahedra.

The distributions of the values of the parameters are given in Figure 10, Figure 11 and Figure 12 and in

Table 7. Mesh quality of the real 3D model after optimization with TetGen program.

Name	Definition	Ideal	Min	Max	Max/Ideal
$SV$	$\sqrt{S}/\sqrt[3]{V}$	$\sqrt[4]{3}\sqrt[6]{72}\approx 2.68433$	$2.68435$	$8.83147$	$3.290$
$ER$	$E_{max}/R$	$2\sqrt{6}\approx 4.8990$	$4.9108$	$605.693$	$123.636$
$EH$	$E_{max}/H_{min}$	$\sqrt{3}/\sqrt{2}\approx 1.2247$	$1.2308$	$301.5059$	$246.127$
$MAXS$	$S_{max}/S$	$1/4=0.25000$	$0.25001$	$0.49804$	$1.992$
$MINS$	$S_{min}/S$	$1/4=0.25000$	$0.00228$	$0.24999$	109.649 ${}^1$

¹ This value is calculated differently—as ideal/min.

.

5. First Comparison of the Models

The quality of all 3D meshes of the models are given in

Table 8. Mesh quality measures of all the models.

Name	Definition	Ideal	Mean	SD	Median	P90	P99	Max
GEOSPHERE:
$SV$	$\sqrt{S}/\sqrt[3]{V}$	$\sqrt[4]{3}\sqrt[6]{72}\approx 2.68433$	$2.84749$	$0.15648$	$2.79918$	$3.03451$	$3.46124$	$4.71468$
$ER$	$E_{max}/R$	$2\sqrt{6}\approx 4.89898$	$7.78874$	$1.82397$	$7.38168$	$9.92162$	$14.42310$	$37.44269$
$EH$	$E_{max}/H_{min}$	$\sqrt{3}/\sqrt{2}\approx 1.22474$	$2.38374$	$0.66151$	$2.24651$	$3.20865$	$4.62071$	$14.62800$
$MAXS$	$S_{max}/S$	$1/4=0.25$	$0.30413$	$0.02705$	$0.30016$	$0.34145$	$0.37957$	$0.44339$
$MINS$	$S_{min}/S$	$1/4=0.25$	$0.19808$	$0.02518$	$0.20173$	0.16303 ${}^1$	0.15112 ${}^1$	0.06385 ${}^1$
GLOBE:
$SV$	$\sqrt{S}/\sqrt[3]{V}$	$\sqrt[4]{3}\sqrt[6]{72}\approx 2.68433$	$2.88477$	$0.20380$	$2.81796$	$3.11881$	$3.74693$	$4.23589$
$ER$	$E_{max}/R$	$2\sqrt{6}\approx 4.89898$	$8.20721$	$2.40916$	$7.56507$	$10.90243$	$18.32399$	$28.12424$
$EH$	$E_{max}/H_{min}$	$\sqrt{3}/\sqrt{2}\approx 1.22474$	$2.54553$	$0.85918$	$2.32799$	$3.59229$	$5.75653$	$11.40720$
$MAXS$	$S_{max}/S$	$1/4=0.25$	$0.30769$	$0.02907$	$0.30364$	$0.34783$	$0.38617$	$0.44492$
$MINS$	$S_{min}/S$	$1/4=0.25$	$0.19443$	$0.02878$	$0.19953$	0.15380 ${}^1$	0.13958 ${}^1$	0.06135 ${}^1$
REAL MODEL:
$SV$	$\sqrt{S}/\sqrt[3]{V}$	$\sqrt[4]{3}\sqrt[6]{72}\approx 2.68433$	$2.950$	$0.329$	$2.821$	$3.347$	$4.284$	$8.831$
$ER$	$E_{max}/R$	$2\sqrt{6}\approx 4.89898$	$9.690$	$6.777$	$7.535$	$15.175$	$37.042$	$605.693$
$EH$	$E_{max}/H_{min}$	$\sqrt{3}/\sqrt{2}\approx 1.22474$	$3.177$	$3.046$	$2.259$	$5.347$	$15.732$	$301.506$
$MAXS$	$S_{max}/S$	$1/4=0.250$	$0.310$	$0.042$	$0.299$	$0.373$	$0.445$	$0.498$
$MINS$	$S_{min}/S$	$1/4=0.250$	$0.188$	$0.048$	$0.205$	0.110 ${}^1$	0.077 ${}^1$	0.002 ${}^1$

¹ In this row the values P90, P99 and max are replaced with P10, P5 and min, respectively.

. The smallest values are always close to the ideal ones; the relative differences was always below 1.23% so those minimal numbers are not presented.

For the simplified models of the bouton (the geosphere and the globe) Table 8 the distance between the mean and the median lay within the range of two standard deviations of the corresponding parameter.

For the real model of the bouton, that takes into consideration the same parameters, the situation was noticeably worse, although it still remained satisfactory. The specified 90th and 99th percentiles were similar for all three models. The maximal values in the case of the real bouton mesh were, however, over one order of magnitude (approximately 30 times) larger than the corresponding values for the simpler models which means that there are some bad quality mesh elements. Nevertheless, their number is extremely low in comparison to the total number of the mesh tetrahedra.

If the simplest structures are to be taken as the reference points for assessing the quality of the real mesh, it must be admitted that the extreme values of the quality parameters (non-standardized) were significantly higher than those reference values. The percentage of the “bad” tetrahedra was, however, low—below 1%, which means that, in general, the meshes were suitable as the input for subsequent optimizations.

6. Further Optimizations Involving Modifications of the Surface and Improvement of Dihedral Angles

The results of examination of the created meshes, particularly the real one, revealed problems with their quality, mainly concerning the distribution of dihedral angles of mesh elements, so additional modifications were needed. That was accomplished by using the Stellar software [29]. Next, the simulation results for the worst mesh were compared with those for the best one.

The Stellar program performs transformations like inserting and deleting vertices, edge contraction, smoothing or removing boundary or near-boundary nodes, removal of boundary edges, multi-face removal and, in addition, compound operations which means combinations of all the above transformations.

Stellar software focuses mainly on correcting dihedral angles of mesh tetrahedra which, according to the author, is crucial for numeric calculations (see [29] and the relevant references therein). The results of improving that parameter are in

Table 9. Mesh quality improvement with Stellar software for all the models (for the real model the worst, “raw”, mesh is also included).

Improv. Stage	Number of Nodes	Number of Faces	Number of Tetrahedra	Ideal DA	Min DA	Max DA	Improv. (Min)	Improv. (Max)
GEOSPHERE:
before	$\mathrm{121,358}$	$\mathrm{63,196}$	$\mathrm{644,999}$	$70.53°$	$5.83°$	$164.92$	—	—
after	$\mathrm{125,915}$	$\mathrm{64,062}$	$\mathrm{472,414}$	$70.53°$	$37.50$	$140.41$	$6.43$	$1.17$
GLOBE:
before	7567	6354	$\mathrm{39,853}$	$70.53°$	$7.90°$	$162.18°$	—	—
after	8850	5096	$\mathrm{26,283}$	$70.53°$	$37.44°$	$139.46°$	$4.74$	$1.16$
REAL MODEL:
worst	$\mathrm{57,060}$	$\mathrm{114,116}$	$\mathrm{194,064}$	$70.53°$	$0.0038°$	$179.915°$	—	—
before	$\mathrm{391,897}$	$\mathrm{296,528}$	$\mathrm{1,996,426}$	$70.53°$	$0.528°$	$175.383°$	$201.06$	$1.03$
after	$\mathrm{572,775}$	$\mathrm{357,578}$	$\mathrm{1,635,827}$	$70.53°$	$19.727°$	$160.273°$	$37.33$	$1.09$

for the three models presented above.

The extreme values of the parameters SV, ER, EH, MAXS and MINS after improving the mesh are given in

Table 10. Mesh quality measures of all the models after improvement.

Name	Definition	Ideal	Mean	SD	Median	P90	P99	Max
GEOSPHERE:
$SV$	$\sqrt{S}/\sqrt[3]{V}$	$\sqrt[4]{3}\sqrt[6]{72}\approx 2.68433$	$2.80129$	$0.06310$	$2.79255$	$2.88640$	$2.97954$	$3.33775$
$ER$	$E_{max}/R$	$2\sqrt{6}\approx 4.89898$	$7.47835$	$1.29978$	$7.26340$	$9.07162$	$11.85719$	$24.99427$
$EH$	$E_{max}/H_{min}$	$\sqrt{3}/\sqrt{2}\approx 1.22474$	$2.31219$	$0.54216$	$2.21765$	$2.98423$	$4.11482$	$9.75960$
$MAXS$	$S_{max}/S$	$1/4=0.25$	$0.30716$	$0.02864$	$0.30297$	$0.34659$	$0.38850$	$0.43438$
$MINS$	$S_{min}/S$	$1/4=0.25$	$0.19643$	$0.02555$	$0.19967$	0.16305 ${}^1$	0.14975 ${}^1$	0.04684 ${}^1$
GLOBE:
$SV$	$\sqrt{S}/\sqrt[3]{V}$	$\sqrt[4]{3}\sqrt[6]{72}\approx 2.68433$	$2.81312$	$0.06964$	$2.80272$	$2.90532$	$3.01900$	$3.21246$
$ER$	$E_{max}/R$	$2\sqrt{6}\approx 4.89898$	$7.76157$	$1.49536$	$7.47481$	$9.66496$	$12.82559$	$19.98282$
$EH$	$E_{max}/H_{min}$	$\sqrt{3}/\sqrt{2}\approx 1.22474$	$2.46626$	$0.66469$	$2.33056$	$3.31735$	$4.71825$	$7.41859$
$MAXS$	$S_{max}/S$	$1/4=0.25$	$0.31461$	$0.03285$	$0.30945$	$0.36080$	$0.40590$	$0.43614$
$MINS$	$S_{min}/S$	$1/4=0.25$	$0.18995$	$0.02885$	$0.19390$	0.15030 ${}^1$	0.13593 ${}^1$	0.06181 ${}^1$
REAL MODEL:
$SV$	$\sqrt{S}/\sqrt[3]{V}$	$\sqrt[4]{3}\sqrt[6]{72}\approx 2.68433$	$2.87951$	$0.14403$	$2.84965$	$3.03650$	$3.37273$	$6.39540$
$ER$	$E_{max}/R$	$2\sqrt{6}\approx 4.89898$	$9.41086$	$6.40490$	$8.34152$	$12.62199$	$24.68540$	$838.46663$
$EH$	$E_{max}/H_{min}$	$\sqrt{3}/\sqrt{2}\approx 1.22474$	$3.17968$	$2.74576$	$2.71390$	$4.63019$	$9.51595$	$383.60024$
$MAXS$	$S_{max}/S$	$1/4=0.250$	$0.32890$	$0.03993$	$0.32307$	$0.38640$	$0.43094$	$0.48228$
$MINS$	$S_{min}/S$	$1/4=0.250$	$0.17597$	$0.03685$	$0.18143$	0.12600 ${}^1$	0.10764 ${}^1$	0.00286 ${}^1$

¹ In this row the values P90, P99 and max are replaced with P10, P5 and min, respectively.

.

7. Simulations

As it has been mentioned above, the last method applied in this paper to compare the quality of meshes was to run example simulations for the different meshes of the same object and comparing the results. The three 3d “real” meshes were chosen: the mesh with no improvement, the one improved by TetGen and the mesh which was additionally optimized with Stellar. The parameters for all three simulations were the same: the size and shape of the surface mesh, the initial density distribution of neurotransmitter (NT), the time step assumed for the calculations and the time interval during which the release of NT took place.

The values determining the results used in this paper were: diffusion coefficient $a=3.0\mathsf{\mu }{\mathrm{m}}^2$/s, permeability coefficient $\alpha =416\mathsf{\mu }{\mathrm{m}}^3$/s, the time step $t=0.0001$ s and the time interval during which the release occurred—from $t=t_0+0.0022$ s to $t=t_0+0.0027$ s.

The criterion used to compare the results was the amount of neurotransmitter released during the first release time interval i.e., from $t=0.0022$ s to $t=0.0027$ s, which embraced 6 consecutive time steps. The additional criterion connected with the first one was the total amount of NT in the bouton.

8. Results

The results are presented in Figure 13. Three plots break down naturally in two groups. The two lower ones, depicting the total amount of neurotransmitter (NT) in the bouton varying in time, are almost identical. The values for the mesh improved with Stellar are only a little lower. This may be attributed to the fact that the total amount of NT was calculated based from the values of the density function in mesh vertices which differed in number and location between meshes. The upper plot (black) refers to the initial unimproved mesh. Its values are remarkably higher and the differences are clearly visible; nearly doubled values at the beginning of the impulse and equal values at the end.

Two out of five values of the real mesh parameters presented in Table 8 increased remarkably after applying Stellar optimization. Namely, the maximal value of ER changed from 605.7 to 838.5 and the maximum of EH—from 301.5 to 383.6. The corresponding 95th percentiles decreased which means that only a very small fraction of the tetrahedra had worse ER or EH quality. The maximal (in case of MINS, minimal) values and the 95th (for MINS, the 5th) percentiles of of all other parameters improved after optimization. The interesting outcome of the optimization is that the only parameter that showed substantial improvement was the one designed for this paper i.e., the SV parameter which changed from 8.83 to 6.40 (a nearly 30% better value). The results confirm that the parameter SV could be calculated instead of six dihedral angles of a tetrahedron. With the improved mesh, simulation of neurotransmitter flow in the examined model should be performed with sufficient precision. The detailed inspection of the accuracy of these calculations will be the subject of further research.

It should be stressed that the mesh analyzed in this paper is quite a complicated one which can be confirmed by the fact that the output mesh quality measured by dihedral angles was only in the range of 19 to 160 degrees. When we take into consideration the statement from [29]: “My implementation usually improves meshes so that no dihedral angle is smaller than ${34}^{°}$ or larger than ${131}^{°}$”, it must be concluded that the mesh considered here is extremely complicated.

9. Effect of Mesh Quality on the Algorithm Performance and Results of Simulations

The algorithm, basically, consists of two modules. The first one is responsible for calculations initialization. Among others, it initializes the vector of the unknown variables i.e., the density vector $\rho$. The second one performs iterative calculations over time and space mesh.

In the second module, the iterative one, the calculations in each iteration are divided into three steps: the first one, during which the equation system is solved without taking release and synthesis into consideration, the second one when the values of the vector $\rho$ are updated with the synthesis, and the third one, that involves the diminishing the value of the neurotransmitter density by the amount resulting from the exocytosis. Within the second said step, solving the equation system must be repeated until the status of each tetrahedron—“producing” or “non-producing”—is established between iterations. This means, among others, that calculations on a single tetrahedron may cause the program get stuck in the loop in case when the tetrahedron has one vertex in the release zone and another one in the synthesis zone. Such vertices exist in the worst, unimproved mesh in the real model.

In general, the comparison of calculations for different meshes is a difficult task. The reason is that vertices of different meshes covering the same object, are usually not identical. Therefore it is not possible to compare results obtained in the same mesh points. In the case of modeling the dynamics of transmitting signals in synapse, however, the situation is convenient because the total release of neurotransmitter into synaptic cleft in the time unit is a direct effect of the run of the processes occurring in the presynaptic bouton. This directness is guaranteed by the fact that the presynaptic bouton, considered as a cybernetic system, has the observability property [6].

Referring to the conducted simulations, the phenomena mentioned above are visible clearly. First of all, for the worst unimproved mesh, the calculations with the synthesis turned on got stuck in a loop because of existence of tetrahedra that intersect both the release area and the synthesis area. Such an example tetrahedron is marked bright lime in Figure 7.

In the case of simulation without synthesis, the total amount of neurotransmitter released into the synaptic cleft differs by a factor of around two, for the worst mesh (Mesh 1) compared to the middle-quality one (Mesh 2) and the best mesh (Mesh 3) which both give practically the same result—see Figure 13 and

Table 11. Release of neurotransmitter during exocytosis phase for three different meshes of the real bouton model [vesicles/time step]. Mesh 1—the worst mesh, Mesh 2—Mesh 1 improved with TetGen, Mesh 3—Mesh 2 improved with Stellar.

Time[s]	Mesh 1	Mesh 2	Mesh 3
0.0022	1.67572458	0.93988859	0.91254579
0.0023	2.52562556	1.25460856	1.25423054
0.0024	1.42968640	0.76737063	0.77887697
0.0025	1.01510141	0.71655252	0.71123427
0.0026	0.80204156	0.59739848	0.60125321
0.0027	0.36672566	0.33349662	0.32721117

. The remarkable difference in the simulated amount of released neurotransmitter means, among others, that even in the case of a simple dynamic model—diffusive-reactive PDE model—the use of a mesh of a poor quality produces results which do not reflect the real solutions. It must be admitted that the shape of the impulse is maintained, and that can be confirmed by comparing impulses conducted via calcium channels measured by other authors [2,3] with the result of our research. In all the papers, including this one, the impulses were V-shaped. It should be also noted that the amount of the released NT, used in this paper, is a very accurate indirect measure of the postsynaptic potential (see Figure 1A in the related work [3]). The peak of the impulse, however, is higher for the worst mesh (Figure 13) which would lead to false conclusions. In genetic research this could even favor the object which produces higher output amplitude as more heterozygous since homozygosity usually means impairing life processes (see a related paper [4], especially Figure 3C therein).

The performed simulations, furthermore, shed—to some extent—light on the problem of the mesh quality in the context of dynamical numerics applications. Namely, in the assessment of the mesh quality the common view is that the dihedral angle between the tetrahedron faces is a significant parameter determining the mesh quality—for the ideal mesh the measure of the dihedral angle between any faces should be $70.53°$. Therefore, the Stellar package is oriented mainly towards improving that parameter The simulations presented in this paper indicate, however, that at least for the diffusive-reactive model the improvement of the parameter connected with the quality of dihedral is not of a great importance—compare the results for Mesh 2 and Mesh 3 in Table 11. Apart from the fact that after applying Stellar package the dihedral angle parameter of the Mesh 3 is significantly better than that of Mesh 2—Table 9—the results for the both meshes are practically the same.

10. Concluding Remarks

As it has been aforementioned in Introduction, the main objective of this paper was to discuss the specifics of the quality measures of the generated mesh with respect to numerical calculations and fidelity of simulations. The obtained results lead to the following conclusions.

• In the situation similar to that considered in this paper, the low quality mesh may alter the results of numerical simulations by a factor greater than two. In the context of modeling transmission of signals in the synapse such a low-fidelity result makes the simulation useless.
• In the contrary to opinions expressed in a few papers, the dihedral angle in the case of our simulations did not prove an important parameter determining the quality of the simulation in the considered case.

The problem of generating the three-dimensional mesh itself is not trivial [30,31]. If the geometry of the simulated object is relatively simple—in our example, a geosphere or a globe with regular internal structures (like the spherically-shaped domain of neurotransmitter supply)—then generation of a 3D mesh of the object does not cause radical problems. The tetrahedral mesh generator like TetGen can easily produce a mesh of a satisfying quality (Table 2, Table 4, and Figure 1). On the other hand, if the geometry of the model is complex, then generation of the mesh is faced with essential difficulties. Even such dedicated software as TetGen is not always capable of dealing with the problem. The quality of the created mesh is significantly worse (see Table 7 as well as Figure 10, Figure 11 and Figure 12). However, even moderately optimized meshes seem to produce satisfying results.