Mathematical Biology - Springer Link

3 downloads 0 Views 321KB Size Report
Oct 9, 2007 - R2C2. U3,. (18). dU3 dt. = 1. R2C3. U2 −. (. 1. R3C3. +. 1. R2C3. ) U3 +. 1. R3C3. E. (19). Above system of ODEs governs the changes of ...
J. Math. Biol. (2008) 56:559–576 DOI 10.1007/s00285-007-0131-5

Mathematical Biology

Model of neurotransmitter fast transport in axon terminal of presynaptic neuron Andrzej Bielecki · Piotr Kalita

Received: 15 May 2007 / Revised: 27 August 2007 / Published online: 9 October 2007 © Springer-Verlag 2007

Abstract In this paper a methodology of mathematical description of the synthesis, storage and release of the neurotransmitter during the fast synaptic transport is presented. The proposed model is based on the initial and boundary value problem for a parabolic nonlinear partial differential equation (PDE). Presented approach enables to express space and time dependences in the process: rate of vesicular replenishment, gradients of vesicular concentration and, through the boundary conditions, the location of docking and release sites. The model should be a good starting point for future numerical simulations since it is based on thoroughly studied parabolic equation. In the article classical and variational formulation of the problem is presented and the unique solution is shown to exist. The model is referred to the model based on ordinary differential equations (ODEs), created by Aristizabal and Glavinovic (AG model). It is shown that, under some assumptions, AG model is a special case of the introduced one. Keywords Synapse · Neurotransmitter · Differential model · Parabolic-type PDE · Reaction diffusion equation Mathematics Subject Classification (2000) 92C20 · 35K57 · 35Q80

A. Bielecki · P. Kalita (B) Institute of Computer Science, Jagiellonian University, Nawojki 11, 30-072 Kraków, Poland e-mail: [email protected] A. Bielecki e-mail: [email protected]

123

560

A. Bielecki, P. Kalita

1 Introduction Modelling of processes in neural system is one of standard topics in studies of biological structures. However, neural systems are such complex structures that, contemporary, models can be only fragmental. In particular the phenomenon of synaptic transmission is intensively investigated both by uncovering the underlying biological mechanisms as well as by mathematical modelling and simulation. The synaptic transport in a chemical synapse is a phenomenon in which important role is played by neurotransmitters—neuroactive substances which are the media of the signal transported between neurons. There are three general approaches to create mathematical models of this phenomenon: – probabilistic models ([2,4], [8, Sect. 7.1.1], [21]), – deterministic models governed by ordinary differential equations, ODEs ([1], [8, Sect. 7.1.4], [11]), – deterministic models governed by partial differential equations, PDEs ([6]-model of axonal transport). In this article the synthesis, storage, diffusion and release of neurotransmitter during fast synaptic transmission is studied. The general aim of this paper is to propose a mathematical model of fast synaptic transmission episode that takes place in the presynapse. The model is based on the nonlinear reaction-diffusion type equation being a parabolic PDE. The model accounts not only for the time variability of considered biological parameters but also their spatial distribution. It should be stressed that in the model of Aristizabal and Glavinovic [1], to which the model presented in this paper is referred, the spatial distribution is taken into account only in simplified way—by means of three pools of vesicles which form a so called compartment model. This paper is an attempt to create the model that more adequately represents the biological reality. In particular: – effective PDE of the model is introduced (Eq. 1), – equation is supported with the initial and boundary conditions (Eqs. 2–5) and the resulting problem is shown to be mathematically correct, i.e. it possesses a unique solution (Theorem 1), – this is shown that the model presented by Aristizabal and Glavinovic [1] is a particular case of the proposed one (Theorem 2). The considered approach belongs to a classical stream of scientific investigations concerning the diffusion motion of particles and vesicles in a neuron. For instance in the paper [6] the authors used parabolic diffusion-like PDE model to describe transport of substances in an axon. In Sect. 2 an initial and boundary value problem for a parabolic PDE that describes the episode of fast synaptic transmission—synthesis, storage and release of synaptic vesicles—is introduced (Eqs. (1)–(5)). First, basing on the neurophysiological foundations (see paragraphs biological foundations, model assumptions and model justification), the classical formulation of the problem is presented. Then the variational formulation is introduced in order to perform the analysis of the model. Initial and boundary conditions and their neurophysiological significance are thoroughly taken

123

Model of neurotransmitter fast transport

561

into consideration. In following Sect. 3 properties of the introduced mathematical model are analyzed. Existence and uniqueness of the solution to the problem is shown. The proof is based on Banach Fixed Point Theorem (compare [5, Chap. 9.2.1, Theorems 1, 2]). In two subsequent sections we refer to the model of the same phenomenon, founded by Aristizabal and Glavinovic (AG model, for abbreviation), presented in details in [1]. First the AG model, basing on a system of ODEs, is briefly recalled. Then, in Sect. 5 it is shown that the AG model can be regarded as a special case of the introduced one described by (1)–(5). Finally concluding remarks are included in last Sect. 6. 2 Mathematical model of fast transmission In this section biological foundations of the modelled phenomenon are presented and discussed. They are followed by mathematical assumptions of the model. Then the formula describing dynamics of a fast transport in a presynaptic neuron is introduced. Biological foundations. During the fast component of synaptic transmission, small molecules (glutamate, acetylcholine, GABA, noradrenaline, serotonin) are released from the docked synaptic vesicles into the synaptic cleft. The molecules are synthesized and packed to vesicles in the axon terminal (terminal bouton) of the presynaptic neuron. The release takes place after the arrival of the brief voltage signal called the action potential to the terminal bouton and inflow of calcium ions through voltage dependent channels. The considered process is depicted in Fig. 1. The presented model describes vesicular storage and release as simple diffusion in the cytoplasm occurring on millisecond time scale and is based on biological evidence and assumptions as given previously ([1,13,14,16–18]).

Fig. 1 Presynaptic episode of the fast transmission. Neurotransmitter is packed to vesicles (1) that diffuse towards the cell membrane (2) where they dock (3). In the response to the inflow of calcium ions (4) triggered by the action potential the neurotransmitter is released from docked vesicles to the synaptic cleft (5)

123

562

A. Bielecki, P. Kalita

Model assumptions. 1. Terminal bouton occupies a fixed 3d domain, the vesicle docking and release sites constitute a fixed part of the domain boundary. 2. The unknown of the model is the concentration of vesicles in the cytoplasm. The unit in which this value is expressed can either be the mass or the quantity of the vesicles or the fraction of cytoplasm volume they occupy. 3. Vesicles diffuse inside the terminal bouton and they are synthesized in some subdomain of the bouton. 4. The efficiency of the vesicle synthesis may either be assumed to be constant or proportional to the difference between the equilibrium concentration (above which the synthesis does not take place) and current concentration. The re-uptake of the released neurotransmitter is not considered in the model. 5. Vesicles do not leave the domain unless the action potential arrives. The arrival of the action potential triggers the possibility of the release of neurotransmitter from vesicles through some fixed period of time. The number of vesicles that release their content in a unit time through the unit area is proportional to the vesicle concentration in the vicinity of the release site. Model justification. In contrast to the most popular approaches in the biological modelling, the model presented in this article takes advantage of the PDEs instead of the ODEs. Such approach, although leads to equations which are more complex and numerical simulations that are more time consuming, allows to take into account the spatial distribution of considered unknown variables (in our case the concentration of vesicles with neurotransmitter as a function of space and time). Such approach allows to enrich the model of the neurotransmission by taking into account the following phenomena: – The vesicles move in the bouton domain so the approach that bases on some abstract pools is only an idealization. Approach presented here can actually capture the movement of vesicles. Concentration of vesicles is assumed to be some continuous function which is justifiable since vesicles can be considered small (their typical diameter is 20–50 nm) with comparison to the bouton (typical diameter of the bouton is several µm). – The ODEs models cannot take into account the fact that the vesicles dock and neurotransmitter is released only at some specific sites on the membrane—that naturally appears in PDE model as ∂Ωd . – Fusion of vesicles to the membrane and resultant exocytosis of neurotransmitter to the cleft causes the decrease of vesicular concentration in the docking sites and this results in the gradients of concentration that trigger the diffusion-type processes that drive vesicles to the membrane. The natural approach to model such processes are diffusion-like PDEs. – In the proposed model synaptic ribbons—specific structures in so called ribbon synapses present for instance in retina cells ([7,12,15]) can be easily taken into consideration. The ribbons form clusters containing 100 or more synaptic vesicles that anchor to the presynaptic membrane close to voltage-gated calcium channels to form an apparently stable depots. This static impression of a depot is confirmed by the observation

123

Model of neurotransmitter fast transport

563

that vesicles diffuse to “hot spots” on the plasma membrane where they are released preferentially ([7,20]). Each cell employs multiple ribbons, ranging from 10 to 100, and consequently, the number of vesicles tethered to all of the ribbons is on the order of 1,000–10,000 ([7,10,15]). In ribbon synapses neurotransmitters are released continuously in contrast to other synapses in which transmitter release is triggered by action potentials. The release of transmitter is then correspondingly short and the neuron signals in an “all or none” fashion. In contrast, neurons that form ribbon synapses do not fire action potentials, but transmit information as graded changes in the membrane potential ([9,12]). It is clear that in ribbons there is significant spatial gradient of vesicular concentration, necessitating the development of models of vesicular storage and release able to account for such spatial gradients. Presented PDE model has such capability. We also remark here that for the case of synaptic ribbons there exists a natural link between the PDE based model and the ODE based compartment model. A compartment, being a pool of vesicles, can then be associated with static spatial domain occupied by the ribbon structure. This idea stands behind the derivation of AG model from the PDE one, as it is presented in Sect. 5. Theoretical verification of the model presented in this paper consists in deriving the AG model used for the same phenomena as the special case of the presented one, that shows that our model is more general. The PDE model in which it is assumed that vesicle pools are static spatial domains is used to derive the AG model by the spatial averaging of the unknowns. It should be stressed that this assumption is a good approximation in the case of synaptic ribbons that are present, e.g. in retina but it could be incorrect in other presynaptic terminals. The model can, however, be developed to take them into account. We remark that the inflow of calcium ions through voltage dependent channels that triggers the release of neurotransmitter after the action potential arrives is not explicitly considered here. The calcium currents, however, are much faster then the movement of vesicles and so in order to take them into account the model would have to operate on two different time scales. In presented model calcium currents are assumed to be immediate and thus they are taken into account only implicitly in the initial and boundary conditions. Model parameters. (i) Ω ⊂ R3 - the domain of the terminal bouton being a sufficiently regular set, n = (n 1 , n 2 , n 3 ) being the outer normal versor to ∂Ω, (ii) Ω3 ⊂ Ω - the domain of neurotransmitter production, (iii) ∂Ωd ⊂ ∂Ω - neurotransmitter release sites on the cell membrane, (iv) β : Ω → R - neurotransmitter source density defined, for example, by β(x) = 0 outside Ω3 and β(x) = βz on Ω3 , (v) ¯ - the balance concentration of vesicles with neurotransmitter in the bouton (new vesicles appear only if the concentration is below this value), (vi) α - the coefficient denoting the rate of neurotransmitter exocytosis, α is the number of vesicles (or molecules) which are released through the unit area of the membrane in unit time by the unit difference of the concentration in the cell and outside the cell (1 action potential activates about 300 vesicles and 1 vesicle contains from 1,000 to 10,000 molecules of neurotransmitter),

123

564

A. Bielecki, P. Kalita

(vii) ai j : Ω → R - the diffusion tensor for the neurotransmitter filled vesicles (assuming that the cytoplasm is isotropic and homogeneous medium, the tensor is diagonal, constant in time and independent on space; the value of all three entries on the diagonal is constant and equal to the diffusion coefficient which, for the acetylcholine, is equal to 300 µm2/sec—see [3]), (viii) τ - the time period through which the neurotransmitter is released from the docked vesicles to the cleft (2-5 µs), (ix) t0 - the arrival moment of the potential. Governing equation. The unknown in the model is the function  : Ω ×[0, T ] → R denoting the concentration of the vesicles with neurotransmitter. The function is the solution to the equation   3  ∂(x, t) ∂ ∂(x, t) ai j (x) + β(x)(¯ − (x, t))+ . = ∂t ∂ xi ∂x j

(1)

i, j=1

Remark. First term on the right hand side of the equation is the classical diffusive term. The second one denotes the neurotransmitter production modelled by the function β. The production is weighed by the term (¯ − (x, t))+ and, in consequence, neurotransmitter is not produced if its concentration is greater then . ¯ This term yields the equation nonlinear. To cope with the nonlinearity it is possible either to remove the weighting term completely or to test its value in some time points and assume that it is constant on the intervals between those points. This approximation will give the equation which is piecewise linear and for small interval lengths it will be good approximation of the nonlinear problem if the concentration does not change fast. The boundary conditions are defined in the following way: – On the whole boundary of Ω apart from the docking sites we assume that no vesicles leave the cytoplasm which leads to the homogeneous Neumann boundary condition 3  i, j=1

ai j

∂(x, t) n i = 0 for (x, t) ∈ (∂Ω\∂Ωd ) × [0, T ]. ∂x j

(2)

– Similarly at the docking sites there is no flow of vesicles apart from the time interval on which the voltage dependent calcium channels are open after the arrival of the action potential 3  i, j=1

ai j

∂(x, t) n i = 0 for (x, t) ∈ ∂Ωd × ([0, t0 ) ∪ (t0 + τ, T ]). ∂x j

(3)

– At the docking sites within the time interval on which ion channels are open (which is triggered by the action potential) the flow of the vesicles outside the

123

Model of neurotransmitter fast transport

565

bouton is proportional to their concentration inside the domain. We assume that vesicle concentration outside Ω is equal to 0. 3  i, j=1

ai j

∂(x, t) n i = −α(x, t) for (x, t) ∈ ∂Ωd × [t0 , t0 + τ ]. ∂x j

(4)

We also assume the initial condition (x, 0) = 0 (x) on Ω.

(5)

Weak formulation. This section presents the derivation of the weak formulation of the defined problem. Such formulation is suitable for the existence, uniqueness and stability analysis as well as for further development of the numerical method used to simulate the modelled process. In order to obtain the weak formulation of the problem let us multiply (1) by the test function v ∈ C ∞ (Ω) and integrate over Ω to get  Ω

   3   ∂ ∂ ∂ ai j v d x + β(¯ − )+ v d x, v dx = ∂t ∂ xi ∂x j i, j=1 Ω

(6)



where the unknown  = (x, t), test function v = v(x) and the parameters ai j = ai j (x), β = β(x). From now on we denote the scalar product in L 2 (Ω) by (·, ·). In the diffusive term we apply the Green formula obtaining (σ is the boundary measure) 

   3  3   ∂ ∂ ∂v ∂ + ai j + ,v = − , ai j n i v dσ +(β ·(−) ¯ , v). (7) ∂t ∂ x j ∂ xi ∂x j i, j=1

i, j=1∂Ω

By the boundary conditions 

   3   ∂ ∂ ∂v ai j − ,v = − , α(t)v dσ + (β · (¯ − )+ , v), ∂t ∂ x j ∂ xi i, j=1

(8)

∂Ωd



where α(t) =

α for t ∈ [t0 , t0 + τ ], 0 otherwise.

(9)

We assume that α > 0 (transmitter flows outwards from the neuron), β ∈ L ∞ (Ω) (efficiency of a neurotransmitter source is bounded), ai j ∈ L ∞ (Ω) and furthermore  (ai j ) satisfies classical coercivity condition i,3 j=1 ai j (x)ξi ξ j ≥ C|ξ |2 for all ξ ∈ R3 almost everywhere in Ω. We also assume that 0 ∈ L 2 (Ω). We denote V = H 1 (Ω) and ·, · denotes the duality pairing between V and V . Finally let W (0, T ) = 2 { ∈ L 2 (0, T ; V ), d dt ∈ L (0, T ; V )}.

123

566

A. Bielecki, P. Kalita

Equation (8) takes the form 

   3   ∂ ∂ ∂v ,v = − , α(t)v dσ + (β · (¯ − )+ , v). ai j − ∂t ∂ x j ∂ xi i, j=1

(10)

∂Ωd

Taking into account that in terminal bouton the concentration of vesicles varies both spatially and temporally in a continuous way, we can give the weak formulation of the considered problem: Find  ∈ W (0, T ) such that (5) and (10) are satisfied for a.e. t ∈ [0, T ] with v ∈ V. (11) 3 Properties of the mathematical model In this section we show the following theorem Theorem 1 The problem (11) has exactly one solution. Proof Denote B[, v] :=

  3   ∂ ∂v ai j + , α(t)v dσ ∂ x j ∂ xi

i, j=1

∂Ωd

and f () := β · (¯ − )+ . Since W (0, T ) is continuously embedded in C([0, T ]; L 2 (Ω)) =: X while proving the existence of a solution Banach Fixed Point Theorem will be used on the space X equipped with the norm

 := max (t) L 2 . 0≤t≤T

Define an operator A : X → X in the following two-step way: Step 1 For  ∈ X define h(t) := f ((t)), t ∈ [0, T ]. It is obvious that f is Lipshitzean. Indeed, since β ∈ L ∞ (Ω) it is obtained | f (1 ) − f (2 )| = |β · (¯ − 1 )+ − β · (¯ − 2 )+ | ≤ β L ∞ · |1 − 2 |. The Lipschitz condition for f implies | f ()| ≤ c · (1 + ||)

123

Model of neurotransmitter fast transport

567

for some positive constant c. It can be easily observed that h ∈ L 2 (0, T ; L 2 (Ω)) T which means that 0 Ω h 2 (t)(x)d xdt < ∞. Indeed T 

T  h (t)(x)d xdt =

0 Ω

T  f ((t))(x)d xdt ≤ 2

2

c2 (1 + 2 (t)(x))d xdt

2

0 Ω

0 Ω

T = c1 + 2c2

1

1

(t) L2 2 dt ≤ c1 + 2c2 T ·  2 < ∞. 0

Step 2 For h ∈ L 2 (0, T ; L 2 (Ω)) we denote by w the solution to the problem 

 dw , v + B[w, v] = (h, v) a.e. in [0, T ], v ∈ V, dt w(0) = 0 .

(12) (13)

To justify the correctness of above definition let us note that the linear problem by which w is defined has exactly one solution which belongs to W (0, T ) (so w ∈ X ). Indeed, B[, v] is coercive (B[v, v] + v 2L 2 ≥ C v 2V for v ∈ V ), continuous (B[, v] ≤ D  V v V for , v ∈ V ) and measurable with respect to the time variable. We can therefore use a classical approach (see for instance [19, Theorem 26.1, p. 397]) basing on Galerkin approximations and energy estimates to find out that w ∈ W (0, T ) is uniquely defined and furthermore the mapping (h, 0 ) → w is continuous from L 2 (0, T ; V ) × L 2 (Ω) to W (0, T ). Define A : X → X by putting w = A[]. Existence of solution. We need to prove that if T > 0 is sufficiently small then A is contractive. Let us then take , ˜ ∈ X and let w = A[], w˜ = A[]. ˜ The mapping w satisfies (12)–(13) for h = f () and w˜ for h˜ = f (). ˜ Subtracting Eq. (12) for w and w˜ from each other we find out 

 d ˜ v) a.e. in [0, T ], v ∈ V, (w − w), ˜ v + B[w − w, ˜ v] = (h − h, dt w(0) − w(0) ˜ = 0.

We take v = w(t) − w(t) ˜ in above equation and use the fact that for v ∈ W (0, T ) d we have v , v = 21 dt

v 2L 2 (see for instance proof of Theorem 25.5, p. 394 in [19]). The coercivity of B yields us d ˜ w − w).

w − w

˜ 2L 2 + 2C w − w

˜ 2V − 2 w − w

˜ 2L 2 ≤ 2(h − h, ˜ dt

(14)

Applying Schwartz inequality to the righthand side of (14) we get d ˜ 2 2. ˜ 2L 2 + h − h

w − w

˜ 2L 2 ≤ 3 w − w

L dt

(15)

123

568

A. Bielecki, P. Kalita

˜ 2 2 and η(t) = w − w

Denoting ξ(t) = h − h

˜ 2L 2 we can write L d η(t) ≤ 3η(t) + ξ(t). dt Applying the Gronwall Lemma (see, e.g. [5, Appendix B.2]) and the fact that η(0) = 0 we obtain ⎛ η(t) ≤ e ⎝η(0) +

t

3t

⎞ ξ(t) dt ⎠ ≤ e3T

0

T

T ξ(t) dt = e

0

f () − f ()

˜ 2L 2 dt.

3T 0

Righthand side of above expression does not depend on t so furthermore T 

w − w

˜ ≤e 2

| f () − f ()| ˜ 2 d xdt.

3T 0 Ω

Now we take advantage of the fact that f is Lipschitzean (with the constant β L ∞ ) T

w − w

˜ ≤e 2

3T

β 2L ∞

 − 

˜ 2L 2 dt ≤ e3T β 2L ∞ T  − 

˜ 2. 0

So we have

A[] − A[]

˜ ≤e

3T 2

1

β L ∞ T 2  −  . ˜

The Lipshitz constant is equal to 0 for T = 0 and is increasing function of T . We can therefore pick T1 > 0 such that this constant is less than 1. By Banach Fixed Point Theorem we get the existence of solution to the Problem (11) on the interval [0, T1 ]. Repeating above technique, starting from point T1 and taking (T1 ) as the initial condition we obtain the solution on the interval [T1 , 2T1 ]. Continuing, after finite number of steps we obtain the solution on the whole [0, T ]. Uniqueness of solution. Let us assume that  and ˜ are two weak solutions of Problem (11). Then we have (15) with w = , h = f () and w˜ = , ˜ h˜ = f () ˜ d

 − 

˜ 2L 2 ≤ 3  − 

˜ 2L 2 + f () − f ()

˜ 2L 2 ≤ (3 + β 2L ∞ )  − 

˜ 2L 2 . (16) dt Again, by the Gronwall Lemma we have  = . ˜

4 AG model of synaptic depression dynamics The model of vesicular storage and release founded by Aristizabal and Glavinovic [1] belongs to a family of compartment models. The authors assume that vesicles are

123

Model of neurotransmitter fast transport

569

located in several pools depending on the proximity to the plasma membrane and the degree of release competence. The dynamics of the vesicular storage is reflected as the transport of the unknown quantity (vesicular density) among three vesicular pools— the immediately available, the small and the large one. Vesicles from immediately available pool (that constitutes about 20% of all vesicles) are released during the action potential and thus the synaptic depression is reflected by depletion of this pool. Immediately available pool is replenished by the flow from the small one which is in turn replenished by the flow from the large pool. The quantities Ci (i = 1, 2, 3) denote the ability of the corresponding pool to store vesicles and the unknown vesicular densities in the pools are denoted by Ui , respectively. Furthermore, Ri1C j , where i, j = 1, 2, 3, denote replenishment rates of various pools and E is the efficiency of the source of vesicles placed inside the large pool. The subscript 1 corresponds to the immediately available pool, whereas 2 and 3 to small and large pool, respectively. Subscript 0 describes the synaptic release channel. When it is closed (during the action potential), the dynamics of the vesicular storage and diffusion is described by the following differential equations   1 1 1 dU1 U1 + =− + U2 , dt R1 C 1 R0 C 1 R1 C 1   dU2 1 1 1 1 U2 + = U1 − + U3 , dt R1 C 2 R2 C 2 R1 C 2 R2 C 2   dU3 1 1 1 1 U3 + = U2 − + E. dt R2 C 3 R3 C 3 R2 C 3 R3 C 3

(17) (18) (19)

Above system of ODEs governs the changes of electric voltage in the circuit consisting of three capacitors and a source connected in parallel—see Fig. 2. Vesicles with neurotransmitter are released during stimulation (i.e. when the switch in the circuit model is closed). The opening of the switch (which corresponds to R0 → ∞) precludes the vesicular release. Then, Eq. (17) becomes 1 dU1 1 =− U1 + U2 . dt R1 C 1 R1 C 1

(20)

If the switch in the loop corresponding to the immediately available pool is closed (which means that the corresponding compartment is turned on) then the circuit is described by Eqs. (17)–(19) whereas otherwise the loop is described by Eqs. (18)–(20). If the Eqs. (17)–(20) are interpreted as a description of the electrical circuit then

Fig. 2 Scheme of electric circuit for the AG model

123

570

A. Bielecki, P. Kalita

Ui denotes potential across the ith capacitor, E is voltage of the source, Ri denotes resistance of the ith loop and Ci is the capacitance of the ith capacitor. 5 The relation between the models Here it is shown that the AG model is a special case of the model introduced in Sects. 2 and 3. AG model is derived by spatial averaging of the unknowns. The space domain is assumed to be the sum of three sets Ω = Ω1 ∪Ω2 ∪Ω3 such that their internals are pairwise disjoint (see Fig. 3), where Ω3 correspond to the reserve pool that contains all the vesicle sources (supp f ⊂ Ω3 ), Ω2 is the “intermediate” readily releasable set through which the vesicles travel towards the cell membrane and Ω1 is the immediately releasable set that contains all the docking and release sites. The boundaries between Ω1 and Ω2 are denoted by ∂Ω21 and between Ω2 and Ω3 by ∂Ω32 respectively and they are assumed to be sufficiently regular. For  being the solution to the presented PDE model we introduce the following notation for the averaged values of the vesicle concentration over the subdomains ¯ 3 (t) =

Ω3

(x, t) d x m(Ω3 )

, ¯ 2 (t) =

Ω2

(x, t) d x m(Ω2 )

, ¯ 1 (t) =

Ω1

(x, t) d x m(Ω1 )

, (21)

where m is the Lebesgue measure. The values ¯ i will be shown to correspond to Ui in the AG model. Let us now summarize the required assumptions: [A] The domain of the presynaptic bouton is assumed to be isotropic so the diffusive term simplifies to the Laplace operator multiplied by the diffusion coefficient a. The domain is assumed to be heterogeneous: the diffusion coefficient depends on the pool (i.e a(x) = ai ⇔ x ∈ Ωi ).

Fig. 3 Domains of the PDE problem

123

Model of neurotransmitter fast transport

571

¯ Such [B] For each t ∈ [0, T ] the solution  to PDE problem belongs to C 2 (Ω). regularity can be obtained under certain assumptions concerning the regularity of the function β in the righthand side of the formula (1) and initial conditions. [C] The concentration of the vesicles never exceeds the value ¯ so the positive part in (1) may be neglected. This is because the model AG which we want to obtain is linear. [D] For almost each point of the boundary ∂Ω32 the averaged values of the solution  are achieved at some points on the normal line. The distances from the boundary to those points are constant in time and, furthermore, they are denoted by l1 (x) and l2 (x). Similar assumption is made for ∂Ω21 and distances are denoted by l3 (x) and l4 (x) respectively (see Fig. 3). This assumption is technical and it is required to approximate the averaged values of the concentration by the difference quotients. [E] The flow through the boundary is proportional not to the concentration in the point of the flow as it is in (4) but is the same in all boundary and proportional to the mean vesicle concentration ¯ 1 . [F] Vesicle source is described by β(x) = βz in Ω3 and β(x) = 0 outside Ω3 . The considered problem is now of the following form ∂(x, t) ∂t (x, t) ∂(x, t) ∂n ∂(x, t) ∂n

= a(x)∆(x, t) + β(x)(¯ − (x, t))+ for (x, t) ∈ Ω × [0, T ], (22) = 0 (x) for (x, t) ∈ Ω × {0},

(23)

= 0 for (x, t) ∈ (∂Ω × [0, T ])\(∂Ωd × [t0 , t0 + τ ]),

(24)

= −α ¯ 1 (t) for (x, t) ∈ ∂Ωd × [t0 , t0 + τ ].

(25)

We show the following theorem. Theorem 2 Under assumptions [A]–[F] if  is the solution to (22)–(25) then ¯ 1 , ¯ 2 , ¯ 3 defined by (21) are the solutions of AG model (17)–(20) for the suitable choice of parameters Ri > 0, Ci > 0, E > 0 and initial conditions Ui (0). Proof Integrating the Eq. (22) over Ω3 the following formula is obtained  Ω3

∂(x, t) d x = a3 ∂t



Ω3

 ∆(x, t) d x +

βz (¯ − (x, t))+ d x.

(26)

Ω3

At the left-hand side it is passed with the integral under the derivative while at the right-hand side the Green Theorem is applied to the first term. By [C] the second term becomes β ¯ z m(Ω3 ) − ¯ 3 (t)βz m(Ω3 ). Thus, m(Ω3 )

d ¯ 3 (t) = a3 dt



∂Ω32

∂(x, t) dσ + β ¯ z m(Ω3 ) − ¯ 3 (t)βz m(Ω3 ). ∂n

(27)

123

572

A. Bielecki, P. Kalita

The normal derivative in (27) is approximated by the difference quotient along the normal line - such approximation is possible owing to [B]. Then we get by [D] a3 d ¯ 3 (t) = dt m(Ω3 )

 ∂Ω32

¯ 2 (t) − ¯ 3 (t) dσ − ¯ 3 (t)βz + β ¯ z. l1 (x) + l2 (x)

(28)

Finally, the following equation is obtained d ¯ 3 (t) a3 = (¯ 2 (t) − ¯ 3 (t)) dt m(Ω3 )

 ∂Ω32

dσ − ¯ 3 (t)βz + β ¯ z. l1 (x) + l2 (x)

(29)

Let us now integrate the Eq. (22) over Ω2  Ω2

∂(x, t) d x = a2 ∂t

 ∆(x, t) d x.

(30)

Ω2

By the Green formula we get m(Ω2 )



d ¯ 2 (t) = a2 dt

∂Ω21



∂(x, t) dσ − a2 ∂n

∂Ω32

∂(x, t) dσ. ∂n

(31)

Now approximating the normal derivatives by the difference quotients along the normal lines we arrive at the ODE   d ¯ 2 (t) dσ dσ a2 a2 = ¯ 1 (t) + ¯ 3 (t) dt m(Ω2 ) l3 (x) + l4 (x) m(Ω2 ) l1 (x) + l2 (x) ∂Ω21



−¯ 2 (t)

a2 ⎜ ⎝ m(Ω2 )

∂Ω32



∂Ω21



dσ + l3 (x) + l4 (x)

∂Ω32



dσ ⎟ ⎠. l1 (x) + l2 (x)

(32)

Finally, integrating the Eq. (22) over Ω1 the following formula is obtained  Ω1

∂(x, t) d x = a1 ∂t

 ∆(x, t) d x.

(33)

Ω1

By the Green formula m(Ω1 )

d ¯ 1 (t) = a1 dt



∂Ω

∂(x, t) dσ − a1 ∂n



∂Ω21

∂(x, t) dσ. ∂n

(34)

Observing the boundary conditions, the first boundary integral in above formula simplifies to the integral on ∂Ωd .

123

Model of neurotransmitter fast transport

573

In the second boundary integral the normal derivative is approximated by the different quotient. Thus a1 d ¯ 1 (t) =− dt m(Ω1 )

 ∂Ωd

a1 α(t)¯ 1 (t) dσ − m(Ω1 ) 

where α(t) =

 ∂Ω21

¯ 1 (t) − ¯ 2 (t) dσ, l4 (x) + l3 (x)

α,

for t ∈ [t0 , t0 + τ ],

0,

for t ∈ [0, t0 ) ∪ (t0 + τ, T ].

(35)

(36)

Finally it is obtained for t ∈ [t0 , t0 + τ ]  d ¯ 1 (t) dσ a1 a1 = −¯ 1 (t) αm(∂Ωd ) − ¯ 1 (t) dt m(Ω1 ) m(Ω1 ) l3 (x) + l4 (x) ∂Ω21  dσ a1 , + ¯ 2 (t) m(Ω1 ) l3 (x) + l4 (x)

(37)

∂Ω21

and, for t ∈ [t0 , t0 + τ ] d ¯ 1 (t) a1 = −¯ 1 (t) dt m(Ω1 )

 ∂Ω21

dσ a1 + ¯ 2 (t) l3 (x) + l4 (x) m(Ω1 )

 ∂Ω21

dσ . l3 (x) + l4 (x) (38)

Summing up, we get the equations (we introduce the symbols l21 := dσ and l32 := ∂Ω32 l1 (x)+l ) 2 (x)



dσ ∂Ω21 l3 (x)+l4 (x)

  d ¯ 1 (t) a1 a1 a1 =− l21 + α(t)m(Ωd ) ¯ 1 (t) + l21 ¯ 2 (t), (39) dt m(Ω1 ) m(Ω1 ) m(Ω1 )   d ¯ 2 (t) a2 a2 a2 a2 = l21 ¯ 1 (t)− l32 + l21 ¯ 2 (t)+ l32 ¯ 3 (t), (40) dt m(Ω2 ) m(Ω2 ) m(Ω2 ) m(Ω2 )   d ¯ 3 (t) a3 a3 = l32 ¯ 2 (t) − l32 + βz ¯ 3 (t) + βz . ¯ (41) dt m(Ω3 ) m(Ω3 ) We observe the correspondence between equations (39)–(41) and (17)–(20): – – – – –

Ui corresponds to ¯ i for i ∈ {1, 2, 3}, i) for i ∈ {1, 2, 3}, respectively, Ci corresponds to m(Ω ai 1 1 , 1 and l32 , respectively, R0 , R1 and R2 correspond to αm(Ω d ) l21 E ¯ in (41), R3 C3 in (19) corresponds to βz  α(t) > 0 if there is action potential [then we retrieve (17)] and 0 if there is no action potential [then we retrieve (20)].

123

574

A. Bielecki, P. Kalita

Fig. 4 The section of the electric circuit (A) and the corresponding ring domain of the heat equation (B)

Initial conditions are obtained by the integration of the initial condition (23) over Ω3 , Ω2 and Ω1 respectively. Interestingly diffusion coefficient corresponds to measure of the subdomain divided by the capacity and geometrical properties of the domain boundary correspond to the resistance. This correspondence becomes clear as we observe Fig. 4. The voltages in the circuit segment presented in A are governed by the equation 1 dUC2 = dt C2



 U R2 U R1 . − R2 R1

On the other hand if we integrate the heat equation ∂ ∂t = a∆ over the ringlike domain Ω2 presented in B we get



a ⎜ d Ω2  d x = ⎝ dt m(Ω2 ) m(Ω2 )

 ∂Ω2

∂ dσ − ∂n



∂Ω1

⎞ ∂ ⎟ dσ ⎠. ∂n

By looking at two equations we see that the diffusion coefficient divided by the measure of the domain corresponds to the inverse of capacitance and the boundary integrals correspond to the currents on the resistors. The capacitor represents the domain Ω2 and the flow of the current through resistors—the flow of the conserved quantity through the boundaries ∂Ω1 and ∂Ω2 . 6 Concluding remarks In this paper the PDE-based model for the process of neurotransmitter filled vesicles diffusion in the terminal bouton of a presynaptic neuron in the presence of action potentials was successfully created. In particular the following aims were achieved: – The proposed PDE model considers both temporal and spatial variability of processes in terminal bouton of presynaptic neuron. – Existence and uniqueness of the solution for proposed model have been proved. – The model is based on the thoroughly studied parabolic equation, which should provide good starting point for future numerical simulations.

123

Model of neurotransmitter fast transport

575

– As the proof of existence and uniqueness is based on the Banach Fixed Point Theorem, numerical simulations can take advantage of iterative mechanism. The convergence of such mechanism is guaranteed. – AG model, by which the parameters of synaptic depression were estimated effectively([1]), turned to be a special case of the model presented in this paper. Considering the introduced model the following continuation of studies can be pointed out: – Investigate the effects of various averagings in the PDE model and their effect on the adequacy of the obtained equations. – Include the re-uptake of neurotransmitter from the cleft back to the presynapse. – Formulate the model for the slow synaptic transmission and conjugate it with the presented model to obtain the system which fully models the kinetics of neuroactive substances in the synapse and the synaptic plasticity. – Run the numerical simulations for the fast transmission, slow one and the conjugate model.

References 1. Aristizabal, F., Glavinovic, M.I.: Simulation and parameter estimation of dynamics of synaptic depression. Biol. Cybern. 90, 3–18 (2004) 2. Boyd, I.A., Martin, A.R.: The end-plate potential in mammalian muscle. J. Physiol. 132, 74–91 (1956) 3. Chang, S., Popov, S.V.: Long-range signaling within growing neuritis mediated by neurotrophin3. Neurobiology 96, 4095–4100 (1999) 4. del Castillo, J., Katz, B.: Quantal components of the end-plate potential. J. Physiol. 124, 560–573 (1954) 5. Evans, L.C.: Partial Differential Equations. American Mathematical Society, Providence (1998) 6. Friedman, A., Craciun, G.: A model of intracellular transport of particles in an axon. J. Math. Biol. 51, 217–246 (2005) 7. von Gersdorff, H., Vardi, E., Matthews, G., Sterling, P.: Evidence that vesicles on the synaptic ribbon of retinal bipolar neurons can be rapidly released. Neuron 16, 1221–1227 (1996) 8. Keener, J., Sneyd, J.: Mathematical Physiology. Springer, New York (1998) 9. Lenzi, D., von Gersdorff, H.: Structure suggests function: the case for synaptic ribbons as exocytotic nanomachines. Bioessays 23, 831–840 (2001) 10. Lenzi, D., Runyeon, J.W., Crum, J., Ellisman, M.H., Roberts, W.M.: Synaptic vesicle populations in saccular hair cells reconstructed by electron tomography. J. Neurosci. 19, 119–132 (1999) 11. Magleby, K.L., Stevens, C.F.: A quantitative description of end-plate currents. J. Physiol. 223, 173–197 (1972) 12. Morgans, C.W.: Neurotransmitter release at ribbon synapses in the retina. Immunol. Cell Biol. 78, 442–446 (2000) 13. Oheim, M., Loerke, D., Stuhmer, W., Chow, R.H.: Multiple stimulation-dependent processes regulate the size of the releasable pool of vesicles. Eur. Biophys. J. 28, 91–101 (1999) 14. Parsons, T.D., Coorsen, J.R., Horstmann, H., Almers, W.: Docked granules, the exocytic burst, and the need for ATP hydrolysis in endocrine cells. Neuron 15, 1085–1096 (1995) 15. Parsons, T.D., Sterling, P.: Synaptic ribbon: conveyor belt or safety belt? Neuron 37, 379–382 (2003) 16. Squire, L.R., Bloom, F.E., McConnel, S.K., Roberts, J.L., Spitzer, L.C., Zigmond, M.J.: Fundamental Neuroscience. Academic, New York (2003) 17. Steyer, J.A., Horstmann, H., Almers, W.: Transport, docking and exocytosis of single granules in live chromaffin cells. Nature 388, 474–478 (1997) 18. Vitale, M.L., Seward, E.P., Trifaró, J.M.: Chromaffin cell cortical actin network dynamics control the size of the release-ready vesicle pool and the initial rate of exocytosis. Neuron 14, 353–363 (1995) 19. Wloka, J.: Partial Differential Equations. Cambridge University Press, New York (1987)

123

576

A. Bielecki, P. Kalita

20. Zenisek, D., Steyer, J.A., Almers, W.: Transport, capture and exocytosis of single synaptic vesicles at active zones. Nature 406, 849–854 (2000) 21. Zucker, R.S., Regehr, W.G.: Short-term synaptic plasticity. Annu. Rev. Physiol. 64, 355–405 (2002)

123