Computer Simulations of Hippocampal Mossy Fiber Cleft Zinc Movements

Abstract

Zinc ions have key regulatory, structural, and catalytic functions and mediate a variety of intra- and intercellular processes. The hippocampal mossy fiber boutons contain large amounts of free or loosely bound vesicular zinc, which can be co-released with glutamate. Zinc can interact with a variety of ionic channels (N-VDCCs, L-VDCCs, KATP), glutamate receptors (AMPA, KA, NMDA 2A, 2B), glutamate transporters (GLAST, EAAT4), and molecules (ATP). The dynamic properties of cleft free, complexed, and total zinc were addressed, considering the known concentration and affinity of various cleft zinc sensitive sites, mainly in the postsynaptic area and in glial cells. The computer model included three different zinc release processes, with short, medium, and long duration, described, like the uptake ones, by alpha functions. The results suggest that, depending on the amount of release, zinc clearance is largely due, either, to zinc binding to NMDA 2A receptor sites or to glial GLAST transporters.

1. Introduction

Zinc is one of the most concentrated trace elements in the brain, being essential for normal cellular function and signaling processes in the central nervous system (CNS) [1, 2, 3]. This system contains very large amounts of chelatable or free zinc [4], mainly in the synaptic vesicles of excitatory nerve terminals [5], essentially in the hippocampal mossy fibers from CA3 area [1, 2]. After release, zinc affects the behavior of several voltage-gated and receptor-operated ionic channels [6, 7, 8, 9, 10, 11, 12]. The action of zinc in different types of receptors and channels depends essentially on two factors: their concentration in the synapses and their affinity for zinc. A clear understanding of the action of zinc in individual binding sites is restricted by the complexity of the synaptic transmission process. To further investigate this zinc role, a computational model was elaborated to describe zinc changes associated with the most important zinc-binding sites in the synaptic cleft between a mossy fiber terminal and a pyramidal cell of the hippocampal CA3 area, as previously reported [12]. Assuming that zinc is co-released with the neurotransmitter glutamate, this model was constructed taking into account previous studies that include computer simulations of glutamate dynamics in the synaptic cleft [13, 14, 15, 16]. In the zinc model [12], the variation of total cleft zinc changes is obtained by subtracting two alpha functions, describing zinc release and zinc uptake. These functions, characterized by rapid climb phases and slower decays, were determined based on the assumed maximum amplitude and rising time values of cleft free zinc concentration. The corresponding parameters were defined taking into account experimental results, from optical and electrophysiological zinc experiments, reporting cleft zinc changes [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27].

The process of zinc clearance from the synaptic cleft may include various actions, such as zinc binding, uptake, and entry into postsynaptic cells. Zinc can bind to a variety of pre- and/or postsynaptic receptors, voltage-dependent ionic channels, glutamate transporters expressed in glial cells, and also free molecules in the cleft medium. One of the most important targets for zinc action is the N-methyl-D-aspartate (NMDA) receptor-binding site with a high affinity for zinc [28]. However, it is present only at low concentrations in the mossy fiber terminals (about 80 nM) [29]. The most abundant zinc-binding site is the GLAST glial glutamate transporter, which is responsible for glutamate removal from the synaptic cleft into the glial cells [30, 31]. Zinc also forms complexes with 2-amino-3-(3-hydroxy-5-methyl-isoxazol-4-yl), propanoic acid (AMPA) and kainic acid (KA) glutamate receptors [29, 32, 33, 34], and potassium-ATP (K_ATP) channels [35, 36], with another type of glutamate transporter, the EAAT4 [37, 38, 39], and with the L- and N-types of voltage-dependent calcium channels (VDCCs) [40, 41, 42]. The concentration, affinity, and kinetics of zinc-binding sites are included in the model and have a very large impact on the behavior of zinc changes. On the other hand, the zinc uptake process is largely unknown, being probably mediated by zinc transporters and/or by zinc movements evoked by the electrochemical gradient [43, 44, 45, 46]. In this model it was assumed that uptake is much slower than release, and a time constant was chosen for the latter process that is much larger than the time constant for the former (see Section 2). With respect to zinc entry into postsynaptic neurons, which may include NMDA receptors, voltage-dependent calcium channels, calcium-permeable AMPA/kainate channels, and the Na/Ca exchanger [17, 18, 32, 34, 47], it is considered that even the strongest stimulation protocol considered in this study is not strong enough to open the postsynaptic zinc permeant routes. In the present model, three different stimulation protocols were considered, named single (zinc release associated with a single stimulus), short, and long (multiple release processes that last for more times). It was considered that those stimulation processes evoked maximum cleft free zinc concentrations of 10 nM, 100 nM, and 1 μM, respectively. These concentrations are close to the values suggested in previous studies performed with similar types of stimulation [11, 19, 21, 27, 48]. Thus, in the present model, only cleft zinc concentrations below or equal to 1 μM were considered; therefore, no zinc enters the postsynaptic region, which may only occur for higher cleft zinc values [18, 20, 45].

The mossy fiber synapses have a very narrow synaptic cleft, measuring less than 20 nm [49, 50, 51]. For that reason, the movement of released glutamate and zinc, with similar free diffusion coefficients, is very rapid, reaching the opposite side of the cleft in a few microseconds. In hippocampal neurons, following an instantaneous release from a vesicle, the decay of glutamate concentration is very fast (tens of microseconds), being reduced to an almost constant value in about 50 μs [16]. Despite the more complex geometry of the mossy fibers [52], it can be assumed that the zinc concentration has a similar time course. For this reason, the diffusion was not included in the model, and it was considered that the cleft free zinc concentration is uniform during the binding process.

2. Methods

2.1 Model equations

The total amount of zinc in the synaptic cleft is, at any one time, partly bound to a number of sites and partly unbound. We shall refer to the concentration of the total amount of zinc in the cleft as [Zn]_T and to the concentration of the unbound or free zinc as [Zn]²⁺. These concentrations may increase due to the release of zinc from the glutamate vesicles when a signal arrives at the presynaptic area and the vesicles open up and pour their contents into the cleft. They may also decrease due to the uptake of zinc into the presynaptic region where it will eventually find its way back into the glutamate vesicles, thus closing the cycle. Therefore, the rate of change of the total concentration of zinc in the cleft is given by

$$\frac{d{\left[\mathit{Zn}\right]}_T}{\mathit{dt}}=R\left(t\right)-U\left(t\right)$$ E1

where $R\left(t\right)$ represents the release, i.e., the rate at which zinc is released from the presynaptic area and enters the cleft, and $U\left(t\right)$ represents the uptake, i.e., the rate at which zinc leaves the cleft and is reabsorbed by the presynaptic region.

As zinc and glutamate are assumed to be released simultaneously, it is reasonable to expect the rates of change of their concentrations to follow the same pattern. Therefore, it is assumed that the function $R\left(t\right)$ will have a fast rising phase, followed by a slow decay, as is known to happen for glutamate. In this study, R(t) and U(t) are described by alpha functions:

$$R\left(t\right)=A_{1}t{e}^{-\frac{t}{{\tau }_1}}$$ E2

where $A_1$ and ${\tau }_1$ are constant values that define the height and position of the peak of the release function, and

$$U\left(t\right)=A_{2}t{e}^{-\frac{t}{{\tau }_2}}$$ E3

where $A_2$ and ${\tau }_2$ are again constant values. We choose ${\tau }_2$ to be much larger than ${\tau }_1$ , as it is well-known that the uptake is usually much slower than the release. These four constants cannot be all independent from each other, as the total concentration of zinc must go back to its resting value so that equilibrium is reached again. $A_1$ , ${\tau }_1$ , and ${\tau }_2$ are chosen to be the independent parameters and $A_2$ to depend on them.

Eq. (1) is easily integrated, yielding

$${\left[\mathit{Zn}\right]}_T={\left[\mathit{Zn}\right]}_{T_r}+A_1{\tau }_1\left[\frac{{\tau }_1}{{\tau }_2}\left(t+{\tau }_2\right)e^{-t/{\tau }_2}-\left(t+{\tau }_1\right)e^{-t/{\tau }_1}\right]$$ E4

where ${\left[\mathit{Zn}\right]}_{T_r}$ denotes the resting value of the total zinc concentration and

$$A_2=A_1{\left(\frac{{\tau }_1}{{\tau }_2}\right)}^2,$$ E5

so as to ensure that ${\left[\mathit{Zn}\right]}_T$ will go back to the basal value ${\left[\mathit{Zn}\right]}_{T_r}$ , as $t\to \infty$ .

Once zinc is released into the cleft, it will react with several different sites, denoted here by $X_i$ , $i=1,\dots ,N$ , thus forming $\mathit{Zn}{X}_i$ complexes. It is assumed that all complexes have 1:1 stoichiometry so that the reaction that takes place can be written as

$$Z{n}^{2+}+X_i\leftrightarrow \mathit{Zn}{X}_i,i=1,\dots ,N.$$ E6

Naturally, the total amount of zinc in the cleft must be equal to the sum of the amount of unbound zinc to the amount of zinc that is bound to all the available sites in the cleft. The total concentration of zinc in the cleft can then be written as

$$\left[\mathit{Zn}\right]{}_T=\left[Z{n}^{2+}\right]+\sum _{i=1}^N\left[\mathit{Zn}{X}_i\right],$$ E7

where $\left[\mathit{Zn}{X}_i\right]$ represents the concentration of the complex $\mathit{Zn}{X}_i$ and the sum extends to all the sites to which zinc can bind.

A similar reasoning can be applied to each binding site $X_i$ . Let ${\left[X_i\right]}_T$ denote the total concentration of $X_i$ in the cleft and $\left[X_i\right]$ the concentration of the free site $X_i$ , that is, the concentration of $X_i$ that is not bound to zinc. Then, assuming that there are no other ions competing with zinc for binding to $X_i$ , one must clearly have

$${\left[X_i\right]}_T=\left[X_i\right]+\left[\mathit{Zn}{X}_i\right],i=1,\dots ,N.$$ E8

The differential equation that describes the dynamics of $\left[\mathit{Zn}{X}_i\right]$ can be obtained from the reaction shown in Eq. (6). Defining $k_{\mathit{on},i}$ and $k_{\mathit{off},i}$ as the association and dissociation rate constants for the reactions involving the site $X_i$ and zinc, this equation becomes

$$\frac{d}{\mathit{dt}}\left[\mathit{Zn}{X}_i\right]=k_{\mathit{on},i}\left[Z{n}^{2+}\right]\left[X_i\right]-k_{\mathit{off},i}\left[\mathit{Zn}{X}_i\right]i=1,\dots ,N.$$ E9

Eq. (8) can be used to remove the variable $\left[X_i\right]$ from Eq. (9) and obtain

$$\frac{d}{\mathit{dt}}\left[\mathit{Zn}{X}_i\right]=k_{\mathit{on},i}\left[{\mathit{Zn}}^{2+}\right]{\left[X_i\right]}_T-\left[\mathit{Zn}{X}_i\right]\left(k_{\mathit{on},i}\left[Z{n}^{2+}\right]+k_{\mathit{off},i}\right)i=1,\dots ,N.$$ E10

Apparently, this equation represents a system of linear, uncoupled, first-order differential equations, which could easily be solved numerically. However, the situation is a little more complicated than that for two reasons.

The first reason is that two parameters required in the differential equation, namely, $k_{\mathit{on},i}$ and $k_{\mathit{off},i}$ , are not known for the reactions of zinc with some of the sites. This issue can be circumvented by assuming that, for those cases, the reactions are so fast, compared to the others, that they quickly adapt to the changes of $\left[Z{n}^{2+}\right]$ . This is done so that the value of $\left[\mathit{Zn}{X}_i\right]$ at any given time is similar to its equilibrium value ${\left[\mathit{Zn}{X}_i\right]}_{\mathit{eq}}$ for that particular free zinc concentration $\left[Z{n}^{2+}\right]$ , which is given by

$$\frac{d}{\mathit{dt}}{\left[\mathit{Zn}{X}_i\right]}_{\mathit{eq}}=0\Rightarrow {\left[\mathit{Zn}{X}_i\right]}_{\mathit{eq}}=\frac{k_{\mathit{on},i}\left[{\mathit{Zn}}^{2+}\right]{\left[X_i\right]}_T}{k_{\mathit{on},i}\left[Z{n}^{2+}\right]+k_{\mathit{off},i}},$$ E11

yielding

$${\left[\mathit{Zn}{X}_i\right]}_{\mathit{eq}}=\frac{\left[Z{n}^{2+}\right]{\left[X_i\right]}_T}{\left[Z{n}^{2+}\right]+k_{D,i}}$$ E12

where the dissociation constant, $k_{D,i}$ , is given by

$$k_{D,i}=\frac{k_{\mathit{off},i}}{k_{\mathit{on},i}}.$$ E13

Thus, if there are known values for $k_{D,i}$ (or IC₅₀ or EC₅₀) for all the reactions in question, the concentrations [ZnX_i ] can be determined using Eq. (12).

The second issue mentioned above is the presence of $\left[Z{n}^{2+}\right]$ in the equations. The concentration of free zinc is an unknown function of time, which will obviously depend on the amount of zinc that is bound to all the sites $X_i$ , so it will depend on the concentrations of all the complexes $\left[\mathit{Zn}{X}_i\right]$ through Eq. (7), which can be written as

$$\left[Z{n}^{2+}\right]=\left[\mathit{Zn}\right]{}_T-\sum _{i=1}^N\left[\mathit{Zn}{X}_i\right].$$ E14

In the model presented here, the total concentration of zinc is a known function of time, so Eq. (14) will couple all the differential Eq. (10) that describe the dynamics of the zinc complexes, as well as the algebraic Eq. (12). This is easy to understand. If more zinc couples to site $X_i$ , there will be less free zinc left to couple to the other sites.

It should be noticed that this model imposes a certain amount of total zinc in the cleft at any specified time. When choosing the parameters of the alpha functions, care must be taken in order not to allow the total amount of zinc in the cleft to be less than the amount of zinc that is bound to the sites X_i at any time.

2.2 Numerical calculations

The task is now to solve the differential Eqs. (10) and the algebraic Eqs. (12), where each one of these equations is coupled to all the others by Eq. (14).

The differential Eqs. (10) can be solved by standard methods, which shall now be described [53]. Let, for the sake of simplicity, the concentrations of the complexes be denoted by f_i(t) and the right-hand side of Eq. (10) be denoted by F_i (t,f₁,…,f_N ). It should be noted that the concentration of free zinc in the right-hand side of that equation will depend on time and on all the complex concentrations, which means that F_i will have to depend also on time and on all the complex concentrations. Eq. (10) will then take the form

$$\frac{d{f}_i}{\mathit{dt}}=F_i\left(t,f_1,\dots ,f_N\right),i=1,\dots ,N.$$ E15

The goal is to integrate this set of equations from t = 0 until t = T for some final time T. In order to achieve that goal, the total time interval T is divided into n equal small time intervals $\mathrm{∆t}$ so that

$$\mathrm{∆t}=\frac{T}{n}.$$ E16

The initial values f_i (0) are assumed to be known. Then, the calculation follows step by step, starting with the evaluation of the values of f_i at time $\mathrm{∆t}$ , next the values at time $2\mathrm{∆t}$ , and so on, until the final time $T=n\mathrm{∆t}$ is reached and the final values are computed. In order to determine the values at time $\left(k+1\right)\mathrm{∆t}$ , where k is an integer, assuming that the values at time $k\mathrm{∆t}$ are known, an expansion of $f_i\left(\left(\mathrm{k}+1\right)\mathrm{∆t}\right)$ in a Taylor series is performed:

$${\mathrm{f}}_{\mathrm{i}}\left(\left(\mathrm{k}+1\right)\mathrm{∆t}\right)={\mathrm{f}}_{\mathrm{i}}\left(\mathrm{k∆t}\right)+{\left[\frac{\mathrm{d}{\mathrm{f}}_{\mathrm{i}}}{\mathrm{dt}}\right]}_{\mathrm{t}=\mathrm{k∆t}}\mathrm{∆t}+\mathrm{O}\left({\mathrm{∆t}}^2\right),$$ E17

where Eq. (15) can be used to yield

$$f_i\left(\left(k+1\right)∆t\right)=f_i\left(k∆t\right)+F_i\left(t,f_1\left(k∆t\right),\dots ,f_N\left(k∆t\right)\right)∆t+O\left({∆t}^2\right).$$ E18

In this equation, $F_i\left(t,f_1,\dots ,f_N\right)$ is evaluated at time $t=k∆t$ with the values of $f_1,\dots ,f_N$ also taken at time $t=k∆t$ : $f_1\left(k∆t\right),\dots ,f_N\left(k∆t\right)$ . Neglecting the term of second order in $∆t$ , Eq. (18) becomes

$$f_i\left(\left(k+1\right)∆t\right)\cong f_i\left(k∆t\right)+F_i\left(t,f_1\left(k∆t\right),\dots ,f_N\left(k∆t\right)\right)∆t.$$ E19

The error in this result is, of course, of second order in $∆t$ .

In order to reach time $t=T=n\mathrm{∆t}$ , it is necessary to take n time steps of length $\mathrm{∆t}$ each. The error in the final result will be the sum of n errors of order ${∆t}^2$ , which makes it of order $O\left(n{∆t}^2\right)=O\left(∆t\right)$ , as n is of order $\frac{1}{∆t}$ , as shown in Eq. (16). Therefore, the error can be made smaller by decreasing the value of $∆t$ , which is equivalent to increasing the value of n. This method is called the Euler method, and the fact that the final error is of first order in $∆t$ makes it a first-order method. It is possible to improve this method by using a more complex calculation at each time step.

In these calculations, the very popular fourth-order Runge–Kutta method was used. In one dimension (only one differential equation), it consists of defining, at time step k

$$g_1=F\left(k∆t,f\left(k∆t\right)\right)∆t$$ E20

$$g_2=F\left(k∆t+\frac{∆t}{2},f\left(k∆t\right)+\frac{g_1}{2}\right)∆t$$ E21

$$g_3=F\left(k∆t+\frac{∆t}{2},f\left(k∆t\right)+\frac{g_2}{2}\right)∆t$$ E22

$$g_4=F\left(k∆t+∆t,f\left(k∆t\right)+g_3\right)∆t.$$ E23

Then, the value of $f\left(\left(\mathrm{k}+1\right)\mathrm{∆t}\right)$ at the next time step will be given by

$$f\left(\left(\mathrm{k}+1\right)\mathrm{∆t}\right)=f\left(\mathrm{k∆t}\right)+\frac{g_1}{6}+\frac{g_2}{3}+\frac{g_3}{3}+\frac{g_4}{6}+O\left({∆t}^5\right).$$ E24

This result can be easily checked by following the lengthy process of expanding all terms in a Taylor series.

After summing all the n increments, the final result will have an error of order ${∆t}^4$ , thus converging much faster than the Euler method. It is easy to extend this result to several dimensions as in the present case.

Usually, the routines that use this method make some sort of quality control of the intermediate results [53]. This can be achieved in several ways, like by comparing a fourth-order result with a fifth-order result. The difference between the two should provide a reasonable estimate for the error. If the error is too small, the routine increases the value of $∆t$ , in order for the calculation to proceed at a faster pace. If the error is too large, then the result is rejected, and a new attempt is made with a smaller step.

This method provides the values of the complex concentrations whose dynamics is described by Eq. (10) after one time step. However, it is also necessary to calculate the values of the other complex concentrations and of the free zinc concentration. This can be achieved by iterating Eqs. (12) and (14) in the following way: using the results obtained by the Runge–Kutta method for the reactions for which the rate constants are known, together with the previous values for the results obtained for the other complex concentrations in Eq. (14), an approximate value for the free zinc concentration is obtained. This value can then be used in Eq. (12), providing new values for the complex concentrations whose rate constants are not known and which are assumed to be in equilibrium. Then this procedure is repeated using these new results in Eq. (14) to obtain again a new value for the free zinc concentration, which is used once more in Eq. (12). The process is repeated until self-consistency is obtained, i.e., until the change of the concentration values after one iteration is negligible.

3. Results

Table 1.

Concentrations of the sites and rate and dissociation constants, IC₅₀ and EC₅₀, for the binding reactions [54, 55]. Source: Reproduced from Quinta-Ferreira et al. [12], with permission from Springer Nature.

Table 2.

Assumed values of free zinc and model parameters for the release (R(t)) and uptake (U(t)) functions, for the single, short, and long zinc release processes. Source: Adapted from Quinta-Ferreira et al. [12], with permission from Springer Nature.

The data indicated in Table 1 allowed a variety of estimates of total, free, and complexed zinc changes at the hippocampal mossy fiber synapses from CA3 area. The studies considered the release processes in a synaptic cleft, both from a single and from multiple vesicles, and also the corresponding uptake. Table 2 shows the data that was used to describe a single vesicle event and two possible multiple vesicle events which are denoted by short and long, where the latter corresponds to a longer time before the maximum concentration of free zinc is obtained, as well as a larger value for that maximum. Table 2 includes also the values of the parameters for the release, R(t), and uptake, U(t), alpha functions defined in Eqs. (2) and (3), for those three cases.

Figure 1 represents the major cellular mechanisms (ionic channels (N-VDCCs, K_ATP and L-VDCCs), glutamate receptors (AMPA, NMDA and KA), and glutamate transporters (GLAST, EAAT4)) forming complexes with cleft free zinc. In this work the high- (NR1a-NR2A) and low (NR1a-NR2B)-affinity NMDA receptor sites will simply be mentioned as NMDA 2A and NMDA 2B, respectively.

Figure 1.

Mossy fiber synaptic components from hippocampal CA3 area. Diagram of the mossy fiber synapse showing synaptic vesicles (SV), glutamate transporters (GLAST and EAAT4), glutamate receptors (NMDA, AMPA, and KA), voltage-dependent calcium channels (N- and L-VDCCs), ATP-sensitive potassium channels (K_ATP), and ATP molecules. The different ions (Zn²⁺, Mg²⁺, Na⁺, and K⁺) are represented by dots and the neurotransmitter glutamate (Glu) by filled triangles. Reproduced from Quinta-Ferreira et al. [12], with permission from Springer Nature.

Since the curves corresponding to the short release process have similar time courses to those of the long release events, as can be seen in Figure 4, only the smaller and faster (single process) and the larger and slower (long process) curves are shown in Figures 2 and 3. Thus, in Figure 2, the upper panels represent superimposed signals of the time derivative of the total zinc, dZn_T/dt; the total zinc, Zn_T; and the free zinc, Zn²⁺, concentrations. Except for the first panel, all the other panels are represented on a semilog scale due to the very large range of the signal amplitudes of the data displayed in Figure 2. This allows for an easier comparison of the initial and maximum concentrations and also of the time courses of the formed complexes. The changes in the dZn_T/dt curves occur very rapidly and are over within 1 and 10 ms, for the single and long processes, respectively. The same curves have maxima at about 100 μs and 1 ms measuring both approximately 55 μM s⁻¹. These curves (see Eqs. (1), (2), and (3)), were obtained subtracting the uptake, U(t), from the release, R(t), alpha functions.

Figure 2.

Modeled zinc changes associated with single and long release processes at hippocampal mossy fiber synapses. (a) The traces represent dZn_T/dt (left), which is given by the difference between the release and the uptake functions, the total (center), and the free (right) zinc concentrations as a function of time. The inset on the left panel shows dZn_T/dt for the single event, in an expanded time scale. (b) Zinc complexes formed with the GLAST (left) glutamate transporter sites and the NMDA 2A (center) and NMDA 2B (right) glutamate receptor sites. In each panel, the smaller and the larger traces are for the single and long release processes, respectively. Note the different scales.

Figure 3.

Zinc complexes for the single and long release processes with time courses similar to those of free zinc. Complexes with unknown on and off-rate constants, assumed to be always at equilibrium with free zinc, formed with (a) AMPA receptors (left), N-VDCCs (center), and K_ATP channels (right). (b) The same type of complexes formed with EAAT4 transporters (left), KA receptors (center), and ATP molecules (right). For both the single (smaller traces) and the long (larger traces) release processes, the signals are displayed by decreasing order of amplitude. Note the different scales.

The lower part of Figure 2 shows again superimposed signals associated with the single and long release processes but for the GLAST, NMDA 2A and NMDA 2B complexes.

In these findings, the peak amplitude of free zinc is about one-half of the peak amplitude of total zinc for the single release process and about two thirds for the long process. Most total and free zinc changes occur within about 5 s, for the single stimulation, and about 15 s, for the long one (Figure 2). The major zinc-binding sites considered in the model are listed in Table 1, where the reaction rate constants and the dissociation constants (or the IC₅₀ or EC₅₀ values), if available from previous works, are also indicated.

The NMDA 2A and NMDA 2B curves, with peaks in the nM range, were built using reaction rate constants. In both cases, the time course for the long process is slower than the time course for the single process (Figure 2).

All the other complexes, with unknown rate constants, are considered to be always in equilibrium with free zinc (see Section 2) and have, thus, identical decreasing shapes to that of free zinc. The dynamics of these complexes, namely, AMPA, N-type VDCCs, K_ATP, EAAT4, KA, and ATP, listed in decreasing order of magnitude, happens in less than 5 s (single process) or 10 s (long process), depending on the intensity of stimulation (Figures 2 and 3). It is interesting to note the faster recovery of the NMDA 2B complex, with respect to NMDA 2A, which is characterized by higher-affinity zinc binding and thus slower dissociation rate, for both the single and long processes.

In the model, based on experimental results, these processes are characterized by an initial free zinc value, $\left[Z{n}^{2+}\right],$ of 1 nM, and by maximum and time-to-peak values of 10 nM and 1 ms for the single stimulation and 1 μM and 100 ms for the long release process. In the single case (and also in the short process), zinc clearance is mainly due to zinc binding to the NMDA 2A and GLAST sites, while in the long stimulation, it is essentially mediated by the formation of GLAST complexes. All the other complexes, characterized by lower affinities, are formed in smaller concentrations, as can be seen in Figures 2 and 3.

An overview of the main results of this study can be observed in Figure 4, for the single, short, and long processes.

Figure 4.

Comparison of the amplitude and time course of the complexes associated with the single, short, and long processes. (a) Superimposed traces of the zinc complexes formed with GLAST and EAAT4 transporters (left), NMDA 2A, NMDA 2B, AMPA, and kainate receptors (center) with K_ATP and N- and L-type VDCCs (right) evoked by a single release event. The latter panel includes also zinc binding to ATP molecules (b) and (c). Similar to (a), but for short (b) and long (c) release processes. Reproduced from Quinta-Ferreira et al. [12], with permission from Springer Nature.

4. Discussion

Mathematical models are a highly valuable tool for the study of synaptic zinc dynamics and in particular of cleft zinc changes. After release, zinc interacts with a variety of pre- and postsynaptic mechanisms which, together with uptake, mediate cleft zinc clearance. In the simpler case, which assumes no postsynaptic zinc entry as considered in this study, all released zinc returns, after some time, to the presynaptic area.

Previous works have suggested that there is no zinc entry in the postsynaptic region for concentrations below 10 μM [18, 20, 27, 56]. As previously reported [12], our model assumes different zinc release events that lead to maximum cleft free zinc concentrations in the range 10 nM to 1 μM, which are below 10 μM. With these values, there should be no zinc entering to the postsynaptic area.

The released zinc can form complexes with various synaptic zinc-binding sites existing mainly on pre- and postsynaptic VDCCs (N- and L-types), K_ATP channels, ionotropic glutamate receptors (AMPA, KA, and NMDA), and also cleft free molecules (ATP) and glial glutamate transporters (EAAT4). The mathematical model considered in this study was designed to obtain estimates of the dynamic behavior of the zinc complexes formed with the mentioned mechanisms [12], assuming previously reported resting and site zinc concentrations and binding constants, which are summarized in Table 1. The formation of the complexes was triggered by three assumed types of stimuli, corresponding to single, short, and long release events, with the maximum and time to peak values, based also on existing experimental findings, indicated in Table 2. In most cases, where the on and off-rate constants were not known, the dissociation constant K_D (or the corresponding IC₅₀ or EC₅₀) values were considered in the equations. In these cases the complexes were considered to attain very rapidly the equilibrium concentration, following each free zinc concentration change. This assumption was made for the following reasons: the steps involved in the formation of a complex include the movement of a zinc ion toward the site, the binding, and the dissociation from it, the latter occurring after a conformational change of the bound target. The value of the on-rate constant is mainly determined by diffusion and binding. The speed of a reaction is largely determined by the diffusion value, considered for all ions as 10⁹–10¹⁰ s⁻¹ [57] and, to some extent, by the rate with which water leaves the ionic solvation sphere. For zinc, this rate and the on-rate constant are considered to be above 10⁷ s⁻¹ and 10⁸ M⁻¹ s⁻¹, respectively [57]. We can now introduce values of on-rate constants between 10⁷ and 10⁸ M⁻¹ s⁻¹, in agreement with the values for the NMDA sites used in the model and the dissociation, IC₅₀ and EC₅₀, values in Table 1. In these cases, the off-rate constants have values in the range 60 s⁻¹ to approximately 7000 s⁻¹, which exceed, by two to four orders of magnitude, the off-rate constant (0.6 s⁻¹) for the NMDA 2A reaction. Thus, since for those complexes unbinding occurs very rapidly, the assumption that they are always in equilibrium with free zinc is justified.

The computational study led to the representation of a variety of time-varying curves to illustrate, for different stimuli, the release minus uptake and the concentrations of total and free zinc, as well as of the various zinc complexes formed. Assuming, in a single release process, a low level of stimulation, causing free zinc to reach 10 nM in the cleft, the predominant complexes are formed with the high-affinity NMDA 2A glutamate receptor sites and with the lower-affinity and highly concentrated GLAST glutamate transporters from glial cells [58]. The concentrations of the other complexes are lower by several orders of magnitude. As mentioned before, all ligands for which the reaction rate constants are not known are assumed to be in equilibrium with free zinc. For this reason all these complexes, and also those associated with the short and long release events, have a similar time course to that of free zinc.

If a 10 times higher free zinc concentration, 100 nM, is attained in the cleft, following the more intense short stimulation, the most abundant zinc complex is now the GLAST one, followed by the NMDA 2A which lasts longer than any of the other complexes. This is due to the much higher affinity of the NMDA 2A sites for zinc. Similar properties apply to the third type of estimated curves, associated with the existence of 1 μM free zinc in the cleft, produced by a longer stimulus. The main differences, with respect to the short stimulus situation, are that the concentration of the NMDA 2A complex remains high for a longer period and all formed complexes have larger amplitudes than in the case of the short stimulus. The fact that externally applied zinc (100 nM) was found to inhibit postsynaptic NMDA currents, at hippocampal CA3 neurons [27], is in agreement with the idea that zinc binds to and inhibits the NMDA receptors. The remaining signals, with much smaller amplitudes, have thus minor or negligible roles in accounting for cleft zinc removal.

The most intense stimulation considered in this study, which assumed a maximum zinc cleft concentration of 1 μM, reveals that, as observed in the short case, GLAST is the complex formed in higher concentration. It is followed by the NMDA 2A complex that saturates when all the corresponding zinc sites (about 40 nM) are bound to zinc. All the other bindings have again much smaller contributions to zinc clearance, especially the complex formed with L-type VDCCs, with a concentration in the fM range.

In conclusion, for a single stimulus, the NMDA 2A high-affinity sites are the most involved in the initial clearance process, while for the stronger stimulations considered (short and long), this role is taken by the highly abundant GLAST complexes. In all cases, uptake has a much slower time course.

At hippocampal CA1 synapses, a single release event may be associated with the release of up to 10 vesicles [58]. In particular, an individual mossy fiber axon has approximately 15 giant boutons (3–5 μm diameter at 150 μm intervals) and 37 release sites (active zones) each [51, 52, 59]. Single boutons contain about 20 active zones at 0.45 μm intervals and 16,000 vesicles, being 1400 ready for release that may be multivesicular, thus leading to very intense release processes [51, 60]. It has been reported that the free or loosely bound zinc and glutamate vesicular concentrations are 1–5 and 60–210 mM, respectively [16, 27, 47]. Let’s assume that zinc and glutamate, which have similar diffusion coefficients [20, 27], are co-released in the same proportion and that, as what happens for glutamate, the cleft released zinc concentration becomes uniform in tens of microseconds [16]. It can be estimated that the volume of a mossy fiber-CA3 cleft (20 nm width, 3 μm radius) is much higher (approximately 17,000 times) than that of a vesicle (40 nm diameter) [27, 49, 50, 51, 61, 62]. The cleft real volume is significantly smaller, by about 25%, if the volume of densely packed conic (20 nm height) dendritic protrusions is subtracted [63]. In this case the discharge of an individual vesicle would lead to the following initial concentrations: 80–400 nM for zinc and 4–20 μM for glutamate. There is a large difference between the radii and, thus, the volumes of the mossy fiber-CA3 and the CA3-CA1 clefts since the latter has a 20 nm width and only 250 nm radius. As a consequence, in the CA1 region, the cleft glutamate range of concentrations is estimated to be around 0.4–2 mM and is thus close to the previously reported ranges of values, 1–5 mM [16] and 0.25–11 mM following an individual vesicle release [64]. Another important issue is the role of the connection between the cleft and extrasynaptic regions with much higher volume [65]. This volume has to be added to that of the cleft, since it forms a large part of the space where the zinc concentration changes occur. A similar fact has been considered for cleft glutamate clearance in neurons of the central nervous system, where the glutamate concentration decreased very rapidly (1–5 ms) 100–500 times [65]. Let us assume again that zinc and glutamate diffuse in a similar way [20, 27]. In this case the concentration range of cleft free zinc changes, evoked by an individual vesicle discharge and after diffusing away from the cleft, will be 0.8–4 nM. If multivesicular release occurs [60], the amount of zinc in the cleft will be significantly larger.

Another estimate of cleft zinc discharge and uptake can be obtained from fluorescent glutamate signals associated with single or repetitive stimulation applied to cultured hippocampal neurons [66]. These authors have found that for the single and short types of stimulation, the maximum concentrations of glutamate were around 0.3 and 0.8 μM, respectively, occurring clearance in less than 1 s, for the single, and 2 s, for the short stimuli. If, as previously estimated, there is about 50 times more glutamate in the cleft, the equivalent maximum concentration range for free zinc will be 6–16 nM.

Previous work has reported that the resting free zinc concentration in the cleft is below 10 nM, meaning that the NMDA 2A sites, which are highly sensitive for zinc, will not become saturated by zinc, and also that the amount of zinc discharged by an individual stimulus does not seem to alter much postsynaptic NMDA currents [27]. Our estimates, for more intense stimulations such as the short and long processes, which are associated with 10 and 100 times more free zinc in the cleft, respectively, than for the single process, indicate that the NMDA 2A sites (K_D = 6 nM) are nearly or fully occupied as reported earlier [10, 27, 28]. For the short and long stimulations, the amount of GLAST complexes is approximately 2 and 10 times more than for NMDA 2A complexes, respectively. Also, only for these stronger stimulation protocols, the much lower-affinity NMDA 2B and AMPA receptors become significantly occupied. The K_ATP channels and the EAAT4 transporters form a reasonable amount of zinc complexes in spite of their smaller affinity, because their concentrations are higher than those of the NMDA sites. As expected, the zinc complexes with lower concentrations and association rates, such as those formed with the kainate receptors and ATP molecules, are only significant following the most intense (long) stimulation. It should also be noticed that the maximum cleft zinc concentration included in the model, 1 μM, does not lead to the inhibition of N- or L-type VDCCs by zinc, since the threshold, half, and almost full blockade concentrations are <5, 69, and 150–200 μM, respectively [40, 42]. For all protocols considered, very small concentrations (in the order of pM-fM) of the VDCC complexes are formed.

The evaluation of the dynamics of synaptic zinc complexes considered in this work contributes to a wider knowledge about synaptic zinc changes. Identifying the main zinc mechanisms involved in mossy fiber zinc clearance is of major importance, considering the potential protective or toxic roles of released zinc at these highly excitable synapses.

Acknowledgments

M.E. Quinta-Ferreira is grateful to Dr. C.C.A.M. Gielen for providing the conditions to start this work. It was funded by the strategic project UID/NEU/04539/2013.