Compartment model of neuropeptide synaptic transport ... - Springer Link

Biol Cybern (2008) 99:443–458 DOI 10.1007/s00422-008-0250-0

ORIGINAL PAPER

Compartment model of neuropeptide synaptic transport with impulse control Andrzej Bielecki · Piotr Kalita · Marian Lewandowski · Marek Skomorowski

Received: 29 October 2007 / Accepted: 25 July 2008 / Published online: 20 September 2008 © Springer-Verlag 2008

Abstract In this paper a mathematical description of a presynaptic episode of slow synaptic neuropeptide transport is proposed. Two interrelated mathematical models, one based on a system of reaction diffusion partial differential equations and another one, a compartment type, based on a system of ordinary differential equations (ODE) are formulated. Processes of inflow, calcium triggered activation, diffusion and release of neuropeptide from large dense core vesicles (LDCV) as well as inflow and diffusion of ionic calcium are represented. The models assume the space constraints on the motion of inactive LDCVs and free diffusion of activated ones and ions of calcium. Numerical simulations for the ODE model are presented as well. Additionally, an electronic circuit, reflecting the functional properties of the mathematically modelled presynaptic slow transport processes, is introduced.

1 Introduction Synaptic transmission in neural systems has usually a chemical character which means that chemical compounds are released to the synaptic cleft from presynaptic boutons. A. Bielecki · P. Kalita (B) · M. Skomorowski Institute of Computer Science, Jagiellonian University, Nawojki 11, 30-072 Kraków, Poland e-mail: [email protected] A. Bielecki e-mail: [email protected] M. Skomorowski e-mail: [email protected] M. Lewandowski Institute of Zoology, Jagiellonian University, Ingardena 6, 30-060 Kraków, Poland e-mail: [email protected]

Neurotransmitter and neuropeptide transport plays a key role in this process. An episode of fast transport, namely the release of neurotransmitter from presynaptic bouton to the synaptic cleft, has been intensively studied, described well on a biological level and modelled using mathematical, cybernetic and electronic models—papers by Aristizabal and Glavinovic (2004) and Kruckenberg and Sandweg (1968) are examples. Mathematical models are usually based on differential equations, both ordinary and partial (Aristizabal and Glavinovic 2004; Friedman and Craciun 2005; Magleby and Stevens 1972) and on difference schemes (Scott and Rusakov 2006). Probabilistic models have been stated as well—see, for instance, Boyd and Martin (1956); del Castillo and Katz (1954); Kalkstein and Magleby 2003 and Keener and Sneyd (1998, Sect. 7.1.1). In recent years the process of slow transport, involving neuropeptides, has been being studied intensively. Fluorescence microscopy allows for both in vitro and in vivo imaging (Han et al. 1999; Levitan et al. 2007; Shakiryanova et al. 2005, 2006) of LDCVs, constituting a large step towards the understanding of underlying mechanisms. It is known that this phenomenon is strictly related to transport of calcium ions and those two processes mutually interact (Greengard 2001): neuropeptide exocytosis and fusion of LDCVs with the membrane are allowed after their activation by ionic calcium. Results concerning a slow excitatory postsynaptic potential (EPSP) and mathematical models of this phenomenon were obtained already in the 1980s, see Adams et al. (1986); Bertrand et al. (2000); however it seems that the kinetics of neuropeptide in a presynaptic bouton has not been so far modelled intensively. This paper is intended to fill this gap. It should be stressed, that the knowledge concerning the slow synaptic transport is, generally, nowadays far incom-

123

444

plete, although it is known that the associated phenomena are crucial for the correct functioning of the neural system. Disorders manifested on the level of slow transmission i.e. neuropeptide synthesis, transport, secretion, hydrolysis as well as in back propagation to the presynaptic cell are the most common causes of such diseases as Alzheimer’s disease (Beal et al. 1986; Gabriel et al. 1996), Parkinson’s disease (Beal et al. 2004), attention-deficit/hyperactivity disorder (Masuo et al. 2004; Oades et al. 2004), schizophrenia (Gabriel et al. 1996; Virgo et al. 1995) and others. Lack of information concerning the mentioned phenomena results, among others, in low efficiency of the medical treatment of the above defects. Theoretical results presented in this article aim to clarify the mechanism of the slow transmission. In the approach described in this paper, two methodologies are used. The first one is based on nonlinear parabolic reaction-diffusion partial differential equations (PDEs), and it follows the approach described in Bielecki and Kalita (2008) for fast transport. The compartment model, built using methodology applied by Aristizabal and Glavinovic (2004), who used a class of linear autonomic ODEs to model neurotransmitter kinetics in a presynaptic bouton, is the second one. The two mentioned methodologies are mutually related. This paper is organized in the following way. In Sect. 2 biological foundations of the phenomena and the formulation of the problem at the modelling level are presented. Then in Sect. 3—the PDE model of considered phenomena is presented. In the succeeding section the ODE model describing the neuropeptide kinetics in a presynaptic bouton is introduced. Section 5 is devoted to the electronic representation of this model. In the next section assumptions concerning the values of the model parameters are discussed. Results of simulations are presented in Sect. 7. Concluding remarks and discussion are given in the last section. Derivation of ODE model from the PDE model as well as that of electronic model, because of their very technical character, are moved to appendices.

2 Problem formulation In this section the neurobiology of a slow neuropeptide transmission, particularly the processes taking place in the presynaptic bouton, which are considered in proposed models, are briefly described (see Greengard (2001); Han et al. (1999); Ng et al. (2003); Shakiryanova et al. (2005, 2006); Levitan (2008) for more detailed description of this phenomenon). Furthermore the motivation for choice of mathematical model is given. Neuropeptides are synthesized in the neuronal cell body where they are packaged in large dense core vesicles (LDCV) and transported to nerve terminals. This is in opposition to the small molecule neurotransmitters. They are stored in the

123

Biol Cybern (2008) 99:443–458

synaptic vesicles which are synthesized inside the bouton and are approximately two times smaller then the LDCVs, which are of size of about 100nm. In the bouton the granules are randomly scattered through the cytoplasm where they slowly accumulate until they fuse with the cell membrane and release their content to synaptic cleft. This release typically occurs slowly (within a time scale of minutes) over long periods in response to repetitive firing. The release of neuropeptides is triggered by Ca2+ influx driven by the large inward driving force through the voltage-gated Ca2+ channels (Scott and Rusakov 2006). Peptidergic neurotransmission is sensitive to patterned electrical activity. It was determined that even the most intense stimuli, e.g. sustained membrane depolarization with Ca2+ , induce slow release of only a fraction of stored neuropeptides, suggesting the presence of a special, releasable pool of secretory vesicles. The nature of this pool is further investigated in other papers (Han et al. 1999; Shakiryanova et al. 2005). Inside the bouton there are immobile (inactive) vesicles and mobile (i.e. rapidly moving) ones. Findings of Han et al. (1999) show that the mobile vesicles diffuse slowly (due to the presence of cytoskeleton) in the bouton and their number decreases by exocytosis. As shown by Shakiryanova et al. (2005) they are scattered randomly in the cytoplasm and their diffusion is undirected. In contrast immobile (inactive) vesicles diffuse much slower and their balance space distribution in the bouton in nonuniform—there are regions in the bouton which are devoid of inactive LDCVs. Mobilization of inactive LDCVs occurs after the repetitive action potential signals—the vesicle activation is triggered by calcium which diffuses to them through the cytoplasm. The process of vesicle diffusion, activation and exocytosis is summarized in the Fig. 1. Partial differential equation models take into consideration the space distributions of unknown variables. Using PDEs to model considered phenomena allows one to take into account the local nature of (possibly anisotropic and heterogeneous) diffusion forces as well as localization of bouton connections with the cell soma and sites of ion channels. Continuous distribution of unknown variables (which is sometimes questionable on microscopic level) and high computational cost of the simulation are two disadvantages of PDE models. Compartment models which use ODEs are based on the assumption that the kinetics of described processes constitute in flows of various substances between different pools. Such flows depend on various factors which are also the part of the model. In our case the flow from the pool of inactive vesicles into the pool of active ones depends on the magnitude of the pool of inactive LDCVs and of Ca2+ (which are the concentrations of inactive LDCVs and of Ca2+ ions in cytoplasm). ODEs reflect only the dependence of modelled quantities (i.e. the concentrations of vesicles in considered pools)

Biol Cybern (2008) 99:443–458

445

Fig. 2 The scheme of an electric circuit for the model of Aristizabal and Glavinovic. Capacitors represent the pools of vesicles, resistors— the diffusion between pools. The source is the synthesis of neurotransmitter, and the switch is responsible for the release

3 PDE model

Fig. 1 After LDCVs filled with neurotransmitter arrive from the cell soma (1) they are captured and immobilized (2). Once the action potential arrives, voltage gated channels open (3) and diffusion and electric forces drive calcium to the cytoplasm, where it activates immobile vesicles (4). Active vesicles diffuse freely through the cytoplasm (5) until they arrive in the vicinity of the membrane (6). There they fuse with the membrane and the neuropeptide is released (7)

over time. Spatial dependencies are neglected or considered only in simplified fashion by associating pools with some space domains. Therefore fewer assumptions are required for ODE models—they are independent on the hypothesis concerning the spacial localization of pools. Furthermore, the most important properties of the process i.e. the activation of vesicles induced by the Ca2+ influx and flows between pools can be still described adequately by the system of ODEs. Summing up, one can say that ODE models are relatively simple and easy to simulate but they focus only on the most crucial aspects of the phenomena, which are the dynamics of the underlying processes. According to the widely accepted compartment methodology used, among others, in Aristizabal and Glavinovic (2004) for the description of the fast transport of neurotransmitters, synaptic depression reflects the depletion of the immediately available pool of transmitter. In modelling kinetics of release, it is often considered that vesicles are localized in several pools depending on the proximity to the plasma membrane and the degree of release competence. Furthermore, only a fraction of all vesicles are available for release. The authors simulated the dynamics of the vesicular storage considering three vesicular pools—the immediately available, the small and the large one. They founded also an electronic circuit model of a secretory cell—see Fig. 2. Vesicles with neurotransmitter are released until the synaptic release channel, represented in the circuit by the switch, is open. If it is closed, what corresponds to the fact that the switch in the circuit model is open, then the neurotransmitter is not transported from the cell.

The model described in this section uses parabolic reaction diffusion PDEs to represent space dependencies (like local diffusion forces driven by concentration gradients or the localization of ion channels and release sites) in the model of LDCVs accumulation, activation, diffusion and release. We make the following assumptions about the problem domains (see also Fig. 3): • The open set Ω ⊂ Rn , n ∈ {2, 3} is the domain of the synaptic bouton. • The open set Ω1
Ω is the part of Ω in which inactive LDCVs accumulate. The set Ω\Ω1 represents the central part of bouton which is devoid of inactive vesicles. • ∂Ω
∂Ω1 i.e. the boundary of the set Ω\Ω1 is contained in the interior of Ω. It means that inactive LDCVs accumulate in the neighborhood of the boundary ∂Ω. • ∂ΩCa
∂Ω is the part of the boundary of Ω in which the calcium channels are located. • ∂Ω N
∂Ω is the part of the boundary of Ω through which inactive LDCVs enter the bouton, ∂Ω N represents the connection of bouton with the soma of the cell. We

Fig. 3 An example of a configuration of sets used in the PDE model. Ω1 is the set in which immobile vesicles accumulate while Ω\Ω1 is devoid of them. The boundary section ∂ΩCa represents the location of calcium channels, while ∂Ω N is the part of the boundary through which new vesicles arrive

123

446

assume that ∂Ω N ∩ ∂ΩCa = ∅ i.e. no calcium channels can be located on the entrance boundary. • ν is the outer normal versor to the boundary of the sets under consideration. We assume that all the sets are sufficiently regular for ν to exist in every point of boundary.

Biol Cybern (2008) 99:443–458

(PN) The problem that describes the dynamics of inactive neuropeptide: n
∂ρ N (x, t) ∂ ij ∂ρ N (x, t) µ N (x) = ∂t ∂x j ∂ xi i, j=1

Thr + −δ(ρCa (x, t) − ρCa ) ρ N (x, t) on Ω1 ,

Let us now list the assumptions used to define the PDEs of the model • The unknowns of the model are ρ N (x, t), ρ A (x, t) and ρCa (x, t), the concentrations of, respectively, inactive LDCVs, active ones and ionic calcium. • Both calcium and active vesicles are allowed to diffuse freely over all bouton domain Ω. The corresponding diffusion coefficients are denoted as µCa and µ A . • The movement of inactive vesicles is allowed only in subdomain Ω1 . Furthermore the diffusion of inactive LDCVs is heterogeneous (i.e. the intensity of diffusion varies between points of Ω1 ) and anisotropic (i.e. some directions of diffusion are more preferred then the others). The difij fusion tensor is denoted as µ N (x), i, j ∈ {1, . . . , n}. • Activation of vesicles takes place only if concentration Thr . For of calcium locally exceeds the threshold value ρCa higher values the intensity of activation depends linearly on calcium concentration. The reaction rate constants for inactive vesicles, active ones and calcium are denoted, respectively, as δ, β and α. • New inactive vesicles are accumulated only if the inactive vesicle concentration locally drops below the value E N . They arrive to the bouton via the boundary part ∂Ω N that defines the connection with axon. The rate at which they arrive is proportional to the difference between E N and their current local concentration. The coefficient of this proportion is given by η. • Calcium is allowed to enter the bouton via the part of its boundary ∂ΩCa that contains the calcium channels. It enters the bouton only if its concentration is below some balance value E Ca and the channels are open after the arrival of action potential signal. The function f (t) defines the capacity of channels. It is equal to 0 if the channels are closed, and if they are open, it is equal to the rate of the flow. • Active vesicles can fuse with the cell membrane anywhere, apart from the connection with the cell soma. The rate at which neuropeptide is exocyted depends linearly on local concentration of LDCVs. This proportion is defined by the constant γ . • Initial distributions of the concentrations ρ N , ρ A and ρCa are known and given, respectively, by ρ N 0 , ρ A0 and ρCa0 . System of three initial and boundary value problems for reaction diffusion PDEs is given as follows:

123

(1)

ρ N (x, t) = 0 on Ω\Ω1 , (2) n
∂ρ N (x, t) ij µ N (x) ν j = η(E N − ρ N (x, t))+ on ∂Ω N , (3) ∂ xi

i, j=1 n
i, j=1

ij

µ N (x)

∂ρ N (x, t) ν j = 0 on ∂Ω1 \∂Ω N , ∂ xi

ρ N (x, 0) = ρ N 0 (x) on Ω1 .

(4) (5)

(PA) The problem that describes the dynamics of active neuropeptide: ∂ρ A (x, t) = µ A ∆ρ A (x, t) ∂t Thr + +β(ρCa (x, t) − ρCa ) ρ N (x, t) on Ω, ∂ρ A (x, t) = −γρ A (x, t) on ∂Ω\∂Ω N , ∂ν ∂ρ A (x, t) = 0 on ∂Ω N , ∂ν ρ A (x, 0) = ρ A0 (x) on Ω.

(6) (7) (8) (9)

(PCa) The problem that describes the calcium dynamics: ∂ρCa (x, t) = µCa ∆ρCa (x, t) ∂t Thr + −α(ρCa (x, t) − ρCa ) ρ N (x, t) on Ω, ∂ρCa (x, t) = f (t)(E Ca − ρCa (x, t))+ on ∂ΩCa , ∂ν ∂ρCa (x, t) = 0 on ∂Ω\∂ΩCa , ∂ν ρCa (x, 0) = ρCa0 (x) on Ω.

(10) (11) (12) (13)

The problems (PCa, PA, PN) are not solved here directly; they are rather used to derive the ODE based compartment model. This approach is similar to that of Bielecki and Kalita (2008) for fast transmission. The unknowns of ODE mode correspond to averaged values of the unknowns of PDE model over the domains of consideration. The following notation for the averaged values of ρ N , ρ A , ρCa over Ω1 and Ω\Ω1 is used.  Ω ρ N (x, t) dx 1 ρ N (t) = 1 , (14) m(Ω1 )  Ω ρ A (x, t) dx 1 ρ A (t) = 1 , (15) m(Ω1 )  Ω\Ω1 ρ A (x, t) dx 2 , (16) ρ A (t) = m(Ω\Ω1 )

Biol Cybern (2008) 99:443–458

 ρCa (t) = 1

ρCa 2 (t) =

Ω1

447

ρCa (x, t) dx

, m(Ω1 )  Ω\Ω1 ρCa (x, t) dx m(Ω\Ω1 )

(17) .

(18)

The defined quantities correspond, respectively, to ρ N , ρ A R , ρ A I , ρCa R and ρCaI in the compartment model presented in the next section.

4 ODE model In this section the model of slow synaptic neuropeptide transmission described by a system of ODEs, based on the compartmental approach, followed by the methodology applied by Aristizabal and Glavinovic (see Aristizabal and Glavinovic 2004), is proposed. Description of variables. In the model we have five pools: one inactive LDCV, two active ones and two calcium. Inactive LDCVs accumulate only in selected regions of the bouton (they are located near the membrane, and the release takes place also in those regions—we denote the corresponding concentrations by subscript R), while both active LDCVs and calcium are allowed to diffuse freely in the whole bouton domain. Let us introduce the following variables: ρ N (t) The concentration of inactive LDCVs with neuropeptide in their accumulation region ρ A R (t) The concentration of activated LDCVs with neuropeptide in the accumulation region of inactive LDCVs ρ A I (t) The concentration of activated LDCVs with neuropeptide in the central region of the bouton devoid of inactive LDCVs ρCa R (t) The concentration of calcium ions in the accumulation region of inactive LDCVs ρCaI (t) The concentration of calcium ions in the central region of the bouton devoid of inactive LDCVs Governing equations. The first equation describes an increase of the vesicles concentration ρ N caused by their transport from the cell soma as well as their activation caused by the inflow of ionic calcium. New vesicles arrive in the bouton only if their concentration inside is lower then E N . The activation yields the decrease of inactive vesicles proportionally to their concentration and the difference between calcium concentration and its threshold value below which the activation does not occur. The equation is dρ N (t) = p1 (E N − ρ N (t))+ dt Thr + −d1 (ρCa R (t) − ρCa ) ρ N (t).

(19)

In the above equation p1 and d1 are constants describing flow rate of vesicles into bouton and the reaction rate constant of the process of vesicle activation. Note that recent findings of Shakiryanova et al. (2006) show that inflow of vesicles is activity dependent. The above equation is based on simplistic assumption that the rate at which new LDCVs come depends only on their concentration. In order to express the activity dependance p1 should depend either explicitly on t or on the capacity of the calcium pool ρCa R . The next equation describes the change of active LDCVs concentration (ρ A R ) in the activation region. This quantity increases due to the activation and decreases due to release into the synaptic cleft. The vesicles are allowed freely to diffuse between the inactive LDCVs accumulation region and central region of the bouton. The direction and the rate of this diffusion depends on the relation between LDCVs concentrations in two regions of the bouton. The equation has the form dρ A R (t) Thr + = d2 (ρCa R (t) − ρCa ) ρ N (t) dt +a1 (ρ A I (t) − ρ A R (t)) − p2 ρ A R (t),

(20)

where d2 is the reaction rate constant for the vesicle activation, a1 represents the diffusion of active LDCVs in the considered pool in the bouton cytoplasm and p2 is the release rate of LDCVs. The following equation describes the change of active vesicles concentration in the central regions of the bouton ρ A I . After the activation this quantity increases because of the diffusion forces that balance the concentration. Later, the central region serves as a reservoir for active vesicles. The equation is dρ A I (t) = −a2 (ρ A I (t) − ρ A R (t)), dt

(21)

where a2 denotes the diffusion constant of active LDCVs in the considered pool. Changes of calcium concentration are taken into account in the last two equations. The first one has three terms that determine the changes of concentration in the inactive LDCVs accumulation region (ρCa R ). The first term represents the inflow of ions through ion channels after the action potential signal (this inflow is proportional to the difference between balance concentration E Ca and current concentration), the second term represents the diffusion of calcium between inactive LDCVs accumulation regions and central regions, while the third term models the use of calcium in the activation process (note that activation takes place above threshold value of calcium concentration and above this value it depends linearly on this concentration above). The equation is

123

448

Biol Cybern (2008) 99:443–458

dρCa R (t) = g(t)(E Ca − ρCa R (t))+ dt +b1 (ρCaI (t) − ρCa R (t)) Thr + ) ρ N (t), −d3 (ρCa R (t) − ρCa

(22)

where b1 denotes diffusion coefficient of calcium in the considered pool, d3 is the reaction rate constant of activation process. The function g is given by the formula  0 during absence of membrane activation, g(t) = p3 > 0 otherwise, where p3 denotes the inflow rate of calcium from the cleft to the bouton when voltage dependent channels are open. The last equation describes the balance of calcium concentration in the central region of the bouton cytoplasm ρCaI . The associated pool fills (by diffusion) with the calcium that is not used to activate immobile LDCVs in pool ρ N , and later it serves as a reservoir for calcium. The equation is dρCaI (t) = −b2 (ρCaI (t) − ρCa R (t)), dt

Table 1 Relation between parameters of PDE model and ODE model. In PDE model the parameters are actual diffusion coefficients, reaction rate parameters and inflow/outflow capacities while parameters of ODE model are lumped ones—i.e. they depend also on domain geometry Parameter type

ODE model

PDE model

Threshold concentration

Thr E Ca , E N , ρCa

Thr E Ca , E N , ρCa

p1

ησ (∂Ω N ) m(Ω1 )

g(t)

µCa f (t)σ (∂ΩCa ) m(Ω1 )

Inflow and outflow

d1

γ µ A σ (∂Ω\∂Ω N ) m(Ω1 ) δ

d2

β

d3

α µ A σ (∂Ω1 \∂Ω) 2hm(Ω1 )

p2

Reaction rate

a1 Diffusion

a2

µ A σ (∂Ω1 \∂Ω) 2hm(Ω\Ω1 )

b1

µCa σ (∂Ω1 \∂Ω) 2hm(Ω1 )

b2

µCa σ (∂Ω1 \∂Ω) 2hm(Ω\Ω1 )

(23)

where b2 denotes the diffusion coefficient of calcium in the considered pool. The structure of mutual dependencies of Eqs. (19-23) is depicted in Fig. 4 and possible distribution of spatial pools is depicted in Fig. 5. Note that, by summing, Eqs. (20) and (21) can be reduced to one equation and, similarly, Eqs. (22) and (23) can also

be reduced to one equation, but such simplification leads to neglecting the diffusion in the considered processes. Relation between the PDE model and the ODE model. To derive the ODE model from the PDE model, we make the following assumption: (H) values of ρ N , ρ A , ρCa in points of Ω1 and Ω\Ω1 can be approximated by the averaged values (14–18).

Fig. 4 The flows between the pools present in the ODE model

Fig. 5 Possible distribution of spatial pools and the placement of processes in the ODE model of slow transport

123

The assumption (H) means that the concentrations do not vary too much in the domains of averaging. This causes that the diffusion forces are considered only between the domains of integration (which are represented by compartments) and not inside them; thus the important feature of PDE model is lost. The space relations, which are encoded explicitly in PDEs, imply certain functional dependencies of the modelled phenomenon (such as for instance the presence of the domain Ω\Ω1 in which the activation does not take place) which are preserved after averaging. Under the assumption (H) the averaged solutions of PDE model are also solutions of ODE model. The proof of this claim is presented in Appendix A. The summary of the relation between constants and parameters of two models (h is the length scale of the domains Ω1 and Ω\Ω1 , m is the Lebesgue measure in Rn and σ is the Lebesgue boundary measure) is presented in the Table 1.

Biol Cybern (2008) 99:443–458

5 Electronic model It has become standard to use electronic circuits to model neural cell functional abilities. Such approach has been intensively exploited to represent processes in both the whole cell and its fragments since early sixties when an electronic model of the whole neuron was proposed by Harmon (1961), see also Tadeusiewicz (1994, Sect. 4). Models of a membrane patch (Nagumo et al. 1964; Keener and Sneyd 1998, Sect. 4.2), a discretized cable model of axonal current (Rall 1969; Keener and Sneyd 1998, Sect. 8.1) and the mentioned AG model can be put as examples. Information processing in the whole neuron and implications for its electronic models were studied by Waxman (1972) and Scott (1995), Sect. 5. It should be also remarked that some studies, concerning possibilities of electronic modelling of signal processing in nervous systems, take into account wide spectrum of neurophysiological phenomena including interneural short-distance and long-distance communication and dendro-dendric coupling. The paper Bielecki et al. (2006b) in which a programmable hybrid (analog-digital and digitalanalog) modules were proposed as an electronic tool for biological neural networks modelling is an example of such approach. Because of the complexity of the system of ODEs (19–23) it is hard to design an electronic circuit that works according to these equations. Therefore a two-step approach according to the top–down methodology was used. In the first step the general structure of the circuit was proposed—see Fig. 6. Electric loops having similar structure to those in AG model were designed in details whereas a fragment of the circuit modelling an LDCVs activation by calcium currents

449

is represented as a black box. This approach corresponds to a classical cybernetic methodology in which a black box is used to describe a system with unknown inner structure but known response for a given input signal (Ashby 1958, chap. 6). In such a way the equations governing the black box behavior were obtained—Eqs. (24–26). Having these equations an electronic circuit governed by them, constituting the inner structure of the black box, is postulated. The second microcontroller models the calcium channel. This is a programmable microcontroller allowing the generation of pulses of any shapes and frequency, in particular those used (E 2 − UC2 )+ . in the Eq. (29), equal to Ch(t) 2 ·R2 The circuit representation of proposed model is depicted in the Fig. 6. The sources E 1 , E 2 , E 3 correspond, respectively, to synthesis of LDCVs, inflow of calcium ions through ion channels and threshold concentration of calcium. Diodes enforce the unidirectional flow of current that corresponds to the fact that in Eqs. (27–29) only positive parts of the terms (E N − ρ N 3 ), (ρCa2 − ρ¯Ca2 ), (E Ca − ρCa1 ) appear. The unknowns ρ N , ρ A I , ρCa R , ρCaI and ρ A R of the model (19–23) which are mean concentrations of LDCVs and Ca2+ ions in the pools are represented by the voltages UC1 , UC2 , UC3 , UC4 and UC5 on the corresponding capacitors in the circuit. Flow resistances between pools (depending on the pools’ geometrical properties) correspond to the resistors. The reaction between ions of calcium and inactive vesicles takes place in a black box which contains electronic microchip that realizes the specific laws that govern input and output voltages and currents. Applying Kirchhoff laws for the presented circuit loops and assuming that the black box is governed by the laws I1 = C1 · const1 · V1 · V2 ,

(24)

I2 = C2 · const2 · V1 · V2 ,

(25)

I3 = C5 · const3 · V1 · V2 ,

(26)

we obtain the following formulae: dUC1 1 = (E 1 − UC1 )+ dt C 1 R1 − const1 (UC2 − E 3 )+ UC1 , dUC5 = const3 (UC2 − E 3 )+ · UC1 dt 1 1 + (UC3 − UC5 ) − UC , C 5 R3 C 5 R5 5

(27)

(28)

1 dUC2 = (UC4 − UC2 ) − const2 (UC2 − E 3 )+ UC1 dt C 2 · R4 Fig. 6 Electric circuit representing the presynaptic episode of slow transmission. The large black box represents the Ca 2+ activation of LDCVs. The laws that govern the box are placed inside. The smaller box represents the pattern of bouton stimulation

+

h(t) (E 2 − UC2 )+ , C 2 · R2

1 dUC3 =− (UC3 − UC5 ) dt C 3 · R3

(29) (30)

123

450

Biol Cybern (2008) 99:443–458 Table 2 Summary of parameters used in simulations. N o means number of LDCVs quantity units, unit volume is an arbitrary unit of cytoplasm volume that can be selected as a unit of reference

Fig. 7 The structure of the black box. A microcontroller is responsible for the multiplication of voltages while the interfaces I, II, III enforce the desired output currents. The interface III is shown in details as an example

and dUC4 1 =− (UC4 − UC2 ). dt R4 · C 4

(31)

Summing up, provided the black box is governed by the laws (24–26), the circuit works according to the Eqs. (27– 31) that exactly correspond to Eqs. (19–23) for the suitable choice of parameters. The black box, in particular, performs multiplication of its input voltages V1 and V2 . This operation can be performed by a programmable microchip. The program should be prepared and implemented on a computer and then the object code should be loaded to the microcontroller. The microcontroller MSC1211 (Texas Instruments Incorporated 2006) can be put as an example. This microchip is equipped, among others, with digital-analog and analog-digital converters as well as the 8051 instructions set compatible core that contains a microprocessor with a memory. Among the inputs - outputs which are used in the presented circuit (see Figs. 6,7) there are analog data and program load inputs as well as digital outputs. The microcontroller communicates with external parts of the circuit via the interfaces I,II,III (see Fig. 7). In particular the interface III, which is shown in details, converts digital output signal into the desired analog current through the integration by the RC loop.

6 The discussion of parameter values The simulations were run for the ODE model (which is simpler and less costly). It should be stressed, however, that the actual values of the ODE model parameters are very hard to determine directly from biological measurements since they encode both physiological parameters and geometric ones coming from the space averaging procedure—see Table 1. For some of the parameters in the model given by Eqs. (19– 23), however, their measured values were taken. Thus E Ca Thr express known concentrations of calcium inside and ρCa and outside the cell (Donahue and Abercrombie 1987). Thre-

123

shold neuropeptide concentration E N is the maximal possible concentration of vesicles in a bouton—which can be expressed in arbitrary units. This calculations are based on the fact that that a total number of LDCVs in single wild-type bouton are about 200 and a diameter of a bouton is about 1.5 µm. Parameters related to initial conditions (see Table 2) refer to the values of the model unknowns at the moment when the biological system observation starts. Those values can vary depending on the synapse state at the moment when the observation starts. Parameters of inflow/outflow, diffusion and reaction (see Table 2) denote the speed of pool replenishment, emptying, and flows between pools. Their values were selected arbitrarily because physical measurements of such quantities are extremely difficult to carry out as functional data concerning the flow of microscopic chemical streams. It seems that nowadays the technology does not allow for measurements of sufficient precision, although it is possible to design such type of experiments. Obtaining such data would allow one to fix the proposed model stronger in physiological reality. Still, the assumed qualitative relations between the parameters levels of magnitudes reflect the known physiological dependencies: diffusion of calcium is much faster then that of neuropeptide, a lot of calcium is required to activate one vesicle, calcium flows into bouton fast, and release of neuropeptide and its replenishment are slow (of which the latter one is slower). The exact values were fit such that the numerical experiment

Biol Cybern (2008) 99:443–458

results are consistent with experimental measurements (see the next section).

7 Simulations for the ODE model Simulations for the ODE model were performed using Matlab, numerical method used to integrate the ODEs was the fourth-order Runge–Kutta scheme with variable time step. Three 500-s simulations of the behavior of the system were run. Initial conditions and values of parameters used in all the simulations are specified in Table 2. In simulations there were applied various patterns of exciting membrane potential, which consisted of 70 Hz bursts of activation as in experiments described in Wong et al. (2008). The main purpose of the first simulation was to verify whether the proposed ODE model is able to express the time changes of LDCVs concentration in Drosophila synaptic bouton measured by Wong et al. (2008) who used a fluorescence method. In their experiment the motor nerve was stimulated by one minute 70 Hz burst of activation. In response rapid release of 20% of neuropeptide was observed. This was followed by a rebound in the following 4 minutes to reach about 90% of the baseline content—see Fig. 1 in the cited paper. In the simulation the stimulation pattern of 70 Hz was represented by function g(t) in (22). After 1 min of continuous stimulation g(t) was set to 0 in order to observe the LDCVs replenishment. Results of simulation are presented in the Fig. 8. Two pools of calcium—ρCa R and ρCa I —almost

Fig. 8 Simulation 1: the bouton was subject to burst of 70 Hz impulses for 60 s after which there were 440 s of inactivity. Within the burst period vesicle pools deplete while slow replenishment is observed in the

451

coincide. The pool Ca R , calcium that flows across the membrane through the ion channel and is used in activation, is changing more rapidly then Ca I which serves as retention basin. The diffusion of calcium is much faster than diffusion of vesicles; therefore Ca I empties almost immediately after Ca R . Because of high values of inflow and reaction coefficients in Eqs. (22) and (23) the 70 Hz pattern reflects almost directly in the changes of calcium concentration. Concurrently a certain smoothing effect can be observed in the Ca I pool. The pool of inactive vesicles N decreases its capacity fast upon stimulation (since the vesicles move to the pool A R ) through the first 60 seconds of experiment and replenishes slowly after the stimulation finishes since the replenishment process is much slower then the activation. The pool asymptotically fills until the balance of concentration is reached. The pool A R increases rapidly within the first 60 seconds corresponding to the strong LDCVs activation. Then the pool empties slowly and continuously by both diffusion to central part of bouton (corresponding to the A I pool) and by the release to the synaptic cleft. These processes proceed until the pool A R is exhausted. Dynamics of the concentration changes in the pool A I is similar to that of A R (but is less rapid since A I plays the role of retention basin). It is confirmed by the observations of Shakiryanova et al. (2005) (see Fig. 3bc in the cited paper) that active vesicles are uniformly distributed in the whole bouton. The total concentration of both active and inactive vesicles ρtotal corresponds to the measured total fluorescence in Wong et al. (2008). It is visible that the results of experiments reflect the observed dynamics of vesicle replenishment—see Fig. 11.

following part of simulation. The concentration ρtotal presented in the last box is the sum ρ N + ρ A R + ρ A I

123

452

Biol Cybern (2008) 99:443–458

Fig. 9 Simulation 2: the bouton was subject to periodic stimulation of 20 s bursts of 70 Hz impulses and 20 s of inactivity. Periodic behavior of all the pools is clearly observed. The concentration ρtotal presented in the last box is the sum ρ N + ρ A R + ρ A I

Fig. 10 Simulation 3: the bouton was subject three randomly chosen bursts of 70 Hz impulses within first 250 s (though which the dominant phenomenon is the release) after which there followed 250 s of inac-

tivity (when the replenishment dominates). The concentration ρtotal presented in the last box is the sum ρ N + ρ A R + ρ A I

In the second simulation (see Fig. 9) the impulse pattern was a repetitive 20 s burst of 70 Hz signals and 20 s silence. The pools Ca R and Ca I reflect the stimulation pattern. Diffusion between them is fast therefore they are close to each other and central pool behaves slightly smoother. The mean concentration of calcium increases within first 200 s of simulation and then the balance is reached between calcium flowing into the bouton via the channels and used on

the activation of vesicles. The inactive vesicles pool depletes in oscillatory manner. Within the burst periods, the release is larger than the refill, and between bursts it is opposite. The general notion is that the dynamics of decrease stabilizes. The pools of active neuropeptide, in turn, increase during bursts and decrease between them. Through first 100 s the increase is larger then the decrease (because the activation is intense then) and, afterwards the process stabilizes with the tendency

123

Biol Cybern (2008) 99:443–458

Fig. 11 Comparison of total quantity of vesicles ρtotal (solid line) obtained in simulation 1 with the experimental data (squares) obtained by Wong et al. (2008) (see Fig. 1), F = 100% in the cited reference corresponds to E N , which was chosen to be equal to 3

of decrease. The total concentration ρtotal decreases since the release is generally faster then replenishment. In the last simulation (Fig. 10) the bouton was exposed to three randomly chosen 70 Hz burst periods within the first 250 s of simulation. Within this period results are similar as in the simulation 2. After the signals stop, within the following 250 s the system slowly returns to its balance, similarly as in the simulation 1.

8 Concluding remarks and comments 1. In the article an effective model of calcium kinetics and release, diffusion, activation and replenishment of neuropeptide-filled LDCVs in the presynaptic bouton controlled by the bursts of electric membrane impulses is presented. Two mathematical models are proposed. The first one is based on PDEs, and it is able to capture the space associated phenomena, namely, the concentration gradient-driven diffusion and location of calcium channels and connection with the cell soma. The mentioned phenomena are modelled by the reaction-diffusion-like system of PDEs similarly as in Bielecki et al. (2006a), Bielecki and Kalita (2008) for the fast transport. The ODE model, in turn, uses the simplistic compartment approach, which neglects the phenomena associated with the spacial dependencies of underlying processes. The compartment model is flexible and more general then the PDE based one because pools can, but do not have to, be associated with space configurations. It is proved in the article that if we associate pools with space configurations then ODE model becomes the special case of PDE one—see Appendix A. Furthermore the simulation results obtained from ODE model comply with the biological measurements—see point 5.

453

2. It should be noticed that in Aristizabal and Glavinovic (2004) the authors remarked that in some cases the sequential arrangement of pools (which they assume) is inadequate. In the models presented in this paper parallel arrangement of pools of Ca and active and inactive neuropeptide was taken into account. Thus the mechanism of activation of neuropeptide by the calcium ions was taken explicitly into account in contrast to the case of fast transport in Aristizabal and Glavinovic (2004). This phenomenon, in contradiction to the case of small vesicles with neurotransmitter, plays the key role in the process and cannot be neglected. 3. The ODE model consists of five pools of active and inactive vesicles and calcium. The model can be reduced to three pools (pools Ca I and A I can be incorporated into Ca R and A R , respectively). Such simplification, however, totally neglects the diffusion and the fact that inactive vesicles accumulate only in some parts of the bouton cytoplasm. 4. Basing on the presented ODE model, the simulations can be run relatively easy. In the article it has been done using the standard software packages (i.e. Matlab/Simulink). 5. The simulation results are physiologically adequate. Some of the simulation results reflect the very recent experimental data, which means that even a simple ODE model is able to capture the actual behavior (see Fig. 11). The model can be further investigated by better fitting the parameter values to the data obtained in experiments under various stimulation routines. 6. The electronic model corresponding to the ODE model has been proposed—see Sect. 5, and its correctness has been shown—see Appendix B. The mathematical models enable studies of biological phenomena on theoretical level. However these phenomena have essential functional aspects. Therefore it is important to investigate if it is possible to implement some functional abilities of biological systems in electronic way in order to use the mentioned abilities in, for instance, robotics. Though the contemporary electronics offer rich possibilities, including digital-analog and analog-digital programming modules it is, at least sometimes, preferable to use the simplest electronic elements—resistors, capacitors, diodes—because of time complexity in programming modules. Therefore we tried to reduce microcontroller usage. 7. The presented model can be further experimentally verified. In particular the changes of calcium concentration can be observed by calcium concentration imaging: a scientific technique usually carried out in research which is designed to show the calcium status of a tissue or medium. By using this method we can observe a rapidly changing intracellular ion concentrations e.g. in presynaptic bouton presented in our model.

123

454

Biol Cybern (2008) 99:443–458

8. The conjugacy of the fast and slow transport, both on the mathematical and electronic levels, is planned to be the continuation of our studies. The extended model which couples the mutual–both inhibitory and excitatory– dependance of slow and fast transport has already been initially elaborated. In the model, which will be the topic of the succeeding article, the synapse that consists of several terminal boutons that interact with each other will be reflected. 9. Note that recent findings of Shakiryanova et al. (2006) show that the bouton acquires vesicles not only from the connection with axon; refill of the vesicle storage in the bouton by activity triggered capture of transiting vesicles was observed. This phenomenon is not considered in the presented PDE and ODE models so far. Taking this phenomenon into account is planned in future work of the authors. Still, the deep problem analysis is required. In particular, the following issues need to be resolved: (a) What are the capture sites? Are they related to the sites of calcium channels? (b) When exactly is the capture triggered? What are the relations between the time period when capture occurs and the time period at which calcium channels are open? (c) How does the capture change with time, is constant throughout the time it occurs? (d) Is the diffusion of acquired vesicles same as the diffusion of the vesicles incoming from the axon? Do they stay immobilized near the membrane or freely diffuse in the cytoplasm?

Appendix A: Derivation of ODE model from the PDE model


n ∂ρ N (x, t) ∂ ij ∂ρ N (x, t) dx = µ (x) dx ∂t ∂x j N ∂ xi i, j=1 Ω1 Ω1  Thr + −δ (ρCa (x, t) − ρCa ) ρ N (x, t) d x.

(32)

Ω1

Since the functions ρ N and ρCa are smooth, exchanging time derivative with the integral, using the Green formula and mean value theorem for integration, we get

∂Ω1

N


∂ρ N (x, t) ij µ N (x) νj ∂ xi i, j=1

Thr + −δm(Ω1 )(ρCa (x1 , t) − ρCa ) ρ N (x1 , t),

123

Thr + −δm(Ω1 )(ρCa (x1 , t)−ρCa ) ρ N (x1 , t).

(34)

Again, using the mean value theorem for integrals we obtain for a point x2 ∈ ∂Ω N , possibly changing with t dρ N 1 (t) = ησ (∂Ω N )(E N − ρ N (x2 , t))+ dt

m(Ω1 )

Thr + −δm(Ω1 )(ρCa (x1 , t)−ρCa ) ρ N (x1 , t).

(35)

Now by assumption (H), we find out that dρ N 1 (t) ησ (∂Ω N ) = (E N − ρ N 1 (t))+ dt m(Ω1 ) Thr + −δ(ρCa 1 (t) − ρCa ) ρ N 1 (t).

(36)

(∂Ω N ) and d1 = δ. Now we This is exactly (19) with p1 = ησm(Ω 1) integrate (6) over Ω1 . We get   ∂ρ A (x, t) d x = µ A ∆ρ A (x, t) d x ∂t

Ω1

 +β

Ω1

Thr + (ρCa (x, t)−ρCa ) ρ N (x, t) d x.

(37)

Ω1

Again ρ A , ρ N and ρCa are sufficiently smooth so we can exchange time derivative with the integral, use the Green formula and mean value theorem for integration. We arrive at  ∂ρ A (x, t) dρ A 1 (t) = µA dσ m(Ω1 ) dt ∂ν Thr + +βm(Ω1 )(ρCa (x3 , t)−ρCa ) ρ N (x3 , t),





∂Ω N

∂Ω1

Let us first integrate (1) over Ω1 . We get

dρ N 1 (t) = m(Ω1 ) dt

where x1 is a point in Ω1 , possibly changing with t and σ is the boundary measure. Using the boundary conditions (3–4) we arrive at  dρ N 1 (t) =η (E N − ρ N (x, t))+ dσ m(Ω1 ) dt

where x3 is a point in Ω1 , possibly changing with t. The boundary ∂Ω1 is the sum of three pairwise disjoint parts ∂Ω1 = ∂Ω N ∪(∂Ω1 \∂Ω)∪(∂Ω\∂Ω N ). Using the boundary conditions (7-8), we arrive at  dρ A 1 (t) ρ A (x, t) dσ m(Ω1 ) = −γ µ A dt  +µ A ∂Ω1 \∂Ω

∂Ω\∂Ω N

∂ρ A (x, t) dσ ∂ν

Thr + +βm(Ω1 )(ρCa (x3 , t)−ρCa ) ρ N (x3 , t).

dσ (33)

(38)

(39)

We apply the mean value theorem for integration to both integrals in above formula, obtaining for some points: x4 ∈

Biol Cybern (2008) 99:443–458

455 A σ (∂Ω1 \∂Ω) . We obtained exactly (21) with a2 = µ2hm(Ω\Ω 1) Similarly we integrate (10) over Ω1 , interchange space integral with time derivative on the left-hand side and use Green formula on the right-hand side to get

∂Ω\∂Ω N and x5 ∈ ∂Ω1 \∂Ω depending on t m(Ω1 )

dρ A 1 (t) = −γ µ A σ (∂Ω\∂Ω N )ρ A (x4 , t) dt

+µ A σ (∂Ω1 \∂Ω)

∂ρ A (x5 , t) ∂ν

Thr + +βm(Ω1 )(ρCa (x3 , t)−ρCa ) ρ N (x3 , t).

(40)

For sufficiently small number h, we find two points x6 ∈ Ω1 and x7 ∈ Ω\Ω1 belonging to the intersection of sphere S(x5 , h) with the line normal to ∂Ω1 \∂Ω and running through x5 . By (H) we take h such that the values ρ A (x6 , t), ρ A (x7 , t) are close to the mean values of ρ A in Ω1 and Ω\Ω1 . Now we approximate the normal derivative in above formula A (x 6 ,t) to get by ρ A (x7 ,t)−ρ 2h m(Ω1 )

dρ A 1 (t) = −γ µ A σ (∂Ω\∂Ω N )ρ A (x4 , t) dt

+µ A σ (∂Ω1 \∂Ω)

 −α

(41)

dρ A dt +

=−

γ µ A σ (∂Ω\∂Ω N ) 1 ρ A (t) m(Ω1 )

Thr + +β(ρCa (t) − ρCa ) ρ N 1 (t). 1

(42)

(∂Ω\∂Ω N ) σ (∂Ω1 \∂Ω) , p2 = γ µ A σm(Ω and Putting a1 = µ A 2hm(Ω 1) 1) d2 = β gives us exactly (20). Let us notice that the flow through the boundary between Ω1 and Ω\Ω1 (which depends on the local concentration gradient on the boundary) is assumed to be proportional to the difference of mean concentrations in two sets. Let us now integrate (6) over Ω\Ω1 . We get   ∂ρ A (x, t) d x = µA ∆ρ A (x, t) d x. ∂t

Ω\Ω1

Ω\Ω1

Similarly as in previous equation, we get m(Ω\Ω1 )

dρ A 2 (t) ∂ρ A (x5 , t) = µ A σ (∂Ω1 \∂Ω) . dt ∂ν

A (x 6 ,t) We use the difference quotient − ρ A (x7 ,t)−ρ to approxi2h ∂ρ A (x5 ,t) mate the normal derivative (note that the direction ∂ν of the outer normal versor is now opposite to the one in the previous equation). Values of ρ A in x6 and x7 are, by (H), approximated by mean values of corresponding domains. We have

dρ A 2 (t) µ A σ (∂Ω1 \∂Ω) =− (ρ A 2 (t) − ρ A 1 (t)). dt 2hm(Ω\Ω1 )

(43)

Furthermore we decompose ∂Ω1 into the sum of three disjoint sets ∂Ω1 = (∂Ω1 \∂Ω) ∪ ∂ΩCa ∪ (∂Ω\∂ΩCa ) and for integrals over last two of them we use boundary conditions (11-12). We have

+

µCa m(Ω1 )

α − m(Ω1 )

µ A σ (∂Ω1 \∂Ω) (ρ A 2 (t) − ρ A 1 (t)) 2hm(Ω1 )

∂Ω1

∂ρCa (x, t) dσ ∂ν

Ω1

Then we use (H) to obtain 1 (t)



Thr + (ρCa (x, t)−ρCa ) ρ N (x, t) d x.

µCa f (t) dρCa 1 (t) = dt m(Ω1 )

ρ A (x7 , t)−ρ A (x6 , t) 2h

Thr + ) ρ N (x3 , t). +βm(Ω1 )(ρCa (x3 , t)−ρCa

dρCa 1 (t) = µCa m(Ω1 ) dt

 ∂Ω1 \∂Ω





(E Ca − ρCa (x, t))+ dσ

∂ΩCa

∂ρCa (x, t) dσ ∂ν

Thr + (ρCa (x, t)−ρCa ) ρ N (x, t) d x.

(44)

Ω1

Now by mean value theorem for integration, there exist points (possibly moving with time) x8 ∈ ∂ΩCa , x9 ∈ ∂Ω1 \∂Ω and x10 ∈ Ω1 such that dρCa 1 (t) µCa f (t)σ (∂ΩCa ) = (E Ca − ρCa (x8 , t))+ dt m(Ω1 ) +

µCa σ (∂Ω1 \∂Ω) ∂ρCa (x9 , t) m(Ω1 ) ∂ν

Thr + −α(ρCa (x10 , t) − ρCa ) ρ N (x10 , t).

(45)

Similarly as in (40) we approximate the normal derivative Ca (x 11 ,t) , where x11 ∈ by the difference quotient ρCa (x12 ,t)−ρ 2h Ω1 , x12 ∈ Ω\Ω1 to obtain dρCa 1 (t) µCa f (t)σ (∂ΩCa ) = (E Ca − ρCa (x8 , t))+ dt m(Ω1 ) +

µCa σ (∂Ω1 \∂Ω) (ρCa (x11 , t) − ρCa (x12 , t)) 2hm(Ω1 )

Thr + ) ρ N (x10 , t). −α(ρCa (x10 , t) − ρCa

(46)

We use (H) to approximate values of unknown functions in specific points of domains of consideration by their corresponding mean values to arrive at

123

456

Biol Cybern (2008) 99:443–458

µCa f (t)σ (∂ΩCa ) dρCa 1 (t) = (E Ca − ρCa 1 (t))+ dt m(Ω1 )

V1 = UC1 it is obtained, from the formula (49) dUC1 1 (E 1 − UC1 )+ − const1 · UC1 · V2 . = dt C 1 · R1

µCa σ (∂Ω1 \∂Ω) (ρCa 2 (t) − ρCa 1 (t)) + 2hm(Ω1 ) Thr + −α(ρCa 1 (t) − ρCa ) ρ N 1 (t).

(47)

σ (∂Ω1 \∂Ω) (t)σ (∂ΩCa ) , b1 = µCa2hm(Ω and d3 = Setting g(t) = µCa f m(Ω 1) 1) α, we have (22). Finally in order to obtain (23), we integrate (10) over Ω\Ω1 , interchange integration with time derivative and use the Green formula, which yields us  ∂ρCa (x, t) dρCa 2 (t) d x = µCa dσ. m(Ω\Ω1 ) dt ∂ν ∂Ω1 \∂Ω

Using mean value theorem for integration and approximating normal derivative by difference quotient, we get dρCa 2 (t) dx m(Ω\Ω1 ) dt = −µCa σ (∂Ω1 \∂Ω)

ρCa (x11 , t) − ρCa (x12 , t) . 2h

(48)

dρCa 2 (t) µCa σ (∂Ω1 \∂Ω) dx = − (ρCa 2 (t) − ρCa 1 (t)). dt 2hm(Ω\Ω1 ) We obtained (23) with b2 =

µCa σ (∂Ω1 \∂Ω) 2hm(Ω\Ω1 ) .

Initial conditions for equations (19-23) are obtained by averaging Eqs. (5, 9, 13) over Ω1 and Ω\Ω1 .

Appendix B: Electronic model justification Consider the loop with the battery E 1 and the capacitor C1 . Because of the diode, the R1 resistor voltage is equal to zero if UC1 > E 1 . Otherwise, by the Kirchhoff voltage law, it is equal to E 1 − UC1 . Thus, introducing the notation  f (x) if f (x) > 0, f + (x) := 0 otherwise, the voltage in the loop is described as (49)

The voltage U R1 = I0 · R1 , where I0 is the R1 resistance current. Applying the Kirchhoff current law we get I0 = I1 + IC1 , where I1 is the current incoming to the black box, dU V1 voltage input - see Fig. 6 – whereas IC1 = C1 · dtC1 is the C1 capacitor current. Using the first black-box law: I1 = C1 · const1 · V1 · V2 and taking into consideration that

123

Considering the loop with the V2 voltage black-box input, battery E 3 and the capacitor C2 , we put V2 = (UC2 − E 3 )+ , into the formula (50) obtaining the formula corresponding to (19) 1 dUC1 = (E 1 − UC1 )+ − const1 (UC2 − E 3 )+ UC1 , dt C 1 R1 (51) where C11R1 corresponds to the constant p in the formula (19) and const1 to d1 . Considering the loop with the capacitors C2 and C4 , we obtain, by the Kirchhoff voltage law UC 2 = UC 4 + I C 4 · R 4 ,

(52) dU

Finally by (H) we have

U R1 = (E 1 − UC1 )+ .

(50)

where the capacitor current IC4 = C4 · dtC4 . This allows us to rewrite the formula (52) in the form explicitly corresponding to the formula (23) 1 dUC4 =− (UC4 − UC2 ), dt R4 · C 4

(53)

where R41·C4 corresponds to the constant b in the formula (23). Let us consider the loop with the capacitor C2 , battery E 2 and resistance R2 . Applying Kirchhoff voltage law we get h(t) · (E 2 − UC2 )+ = I˜2 · R2 ,

(54)

where I˜2 is R2 resistance current and h(t) is a pulse function programmed by the microcontroller. Because, on one hand, the black-box input current I2 is given, by Kirchhoff current law, as I2 = I˜2 −IC2 −IC4 and on the other hand, by the blackbox law, I2 = C2 · const2 · V1 · V2 and recalling that V1 · V2 = (UC2 − E 3 )+ · UC1 , we obtain the formula corresponding to the formula (22) 1 dUC2 = (UC4 − UC2 ) − const2 (UC2 − E 3 )+ UC1 dt C 2 · R4 h(t) + (E 2 − UC2 )+ , (55) C 2 · R2 where

h(t) C2 ·R2

corresponds to g(t) in the formula (22) whereas corresponds to b and const2 to d3 . Considering the voltage drop in the C5 − C3 capacitors loop, we obtain UC5 − UC3 − IC3 R3 = 0 and, as a consequence, the formula corresponding to (21) 1 C2 ·R4

1 dUC3 =− (UC3 − UC5 ). dt C 3 · R3

(56)

Biol Cybern (2008) 99:443–458

457

Applying the Kirchhoff current law to the loop attached to the black-box output, we obtain I 5 = I 3 − IC 3 − IC 5 = C5 · const3 · (UC2 − E 3 )+ · UC1 −C3

dUC3 dUC5 −C5 , dt dt

where I5 is R5 resistor current and IC3 , IC5 are C3 and C5 capacitor currents respectively. Considering the loop with the capacitor C5 and resistor R5 , we have UC5 = I5 R5 and, as a consequence UC5 = R5 C5 · const3 (UC2 − E 3 )+ · UC1 − C3 R5 −C5 R5

dUC5 . dt

dUC3 dt (57)

Substituting (56) to (57), we obtain formula corresponding to (20) dUC5 1 = const3 (UC2 − E 3 )+ · UC1 + (UC3 − UC5 ) dt C 5 R3 1 − UC . (58) C 5 R5 5 Acknowledgments The authors would like to thank the anonymous reviewers for the valuable comments that greatly helped to improve the article.

References Adams PR, Jones SW, Pennefather P, Brown DA, Koch C, Lancaster B (1986) Slow synaptic transmission in frog sympathetic ganglia. J Exp Biol 124:259–285 Ashby WR (1958) An Introduction to Cybernetics. Chapman & Hall Ltd, London Aristizabal F, Glavinovic MI (2004) Simulation and parameter estimation of dynamics of synaptic depression. Biol Cybern 90:3–18 Beal MF, Mazurek MF, Chattha GK, Svendsen CN, Bird ED, Martin JB (1986) Neuropeptide Y immunocreativity is reduced in cerebral cortex in Alzheimer’s disease. Ann Neurol 20:282–288 Beal MF, Clevens RA, Mazurek MF (2004) Somatostatin and neuropeptide Y immunoreactivity in Parkinson’s disease dementia with Alzheimer’s changes. Synapse 2(4):463–467 Bertrand PP, Thomas EA, Kunze WAA, Bornstein JC (2000) A simple mathematical model of second-messenger mediated slow excitatory postsynaptic potentials. J Comput Neurosci 8:127–142 Bielecki A, Kalita P, Lewandowski M (2006a) Model of fast synaptic transmission. In: Proceedings of the 12th national conference on application of mathematics in biology and medicine. Kraków, Poland, pp 13–18 Bielecki A, Skomorowski M, Wo´zniak R (2006b) Electronic neural networks modelling. In: Cader A, Rutkowski L, Tadeusiewicz R, ˙ Zurada J (eds) Artificial intelligence and soft computing. Series: challenging problems of science—computer science, Bolc L (ed) Academic Publishing House EXIT, Warszawa, Poland, pp 8–14 Bielecki A, Kalita P (2008) Model of neurotransmitter fast transport in axon terminal of presynaptic neuron. J Math Biol 56(4):559–576 Boyd IA, Martin AR (1956) The end-plate potential in mammalian muscle. J Physiol 132:74–91

del Castillo J, Katz B (1954) Quantal components of the end-plate potential. J Physiol 124:560–573 Donahue BS, Abercrombie RF (1987) Free diffusion coefficient of ionic calcium in cytoplasm. Cell Calcium 8:437–448 Friedman A, Craciun G (2005) A model of intracellular transport of particles in an axon. J Math Biol 51:217–246 Gabriel SM, Davidson M, Haroutunian V, Powchik P, Bierer LM (1996) Neuropeptide deficits in schizophrenia vs. Alzheimer’s disease cerebral cortex. Biol Psychiatry 39:82–91 Greengard P (2001) The neurobiology of slow synaptic transmission. Science 294:1024–1030 Han W, Ng YK, Axelrod D, Levitan ES (1999) Neuropeptide release by efficient recruitment of diffusing cytoplasmic secretory vesicles. Proc Natl Acad Sci USA 96:14577–14582 Harmon LD (1961) Studies with artificial neuron—properties and functions of an artificial neuron. Biol Cybern 1:89–101 Kalkstein JM, Magleby KL (2003) Augmentation increases vesicular release probability in the presence of masking depression at the frog neuromuscular junction. J Neurosci 24:4127–4134 Keener J, Sneyd J (1998) Mathematical physiology. Springer, New York Kruckenberg P, Sandweg R (1968) An analog model for acetylocholine release by motor nerve endings. J Theoret Biol 19:327–332 Levitan ES, Lanni F, Shakiryanova D (2007) In vivo imaging of vesicle motion and release at the Drosophila neuromuscular junction. Nat Protocols 2:1117–1125 Levitan ES (2008) Signaling for vesicle mobilization and synaptic plasticity. Mol Neurobiol 37(1):39–43 Magleby KL, Stevens CF (1972) A quantitative description of end-plate currents. J Physiol 223:173–197 Masuo Y, Morita M, Oka S, Ishido M (2004) Motor hyperactivity caused by a deficit in dopaminergic neurons and the effects of endocrine disruptors: a study inspired by the physiological roles of PACAP in the brain. Regul Pept 123:225–234 Nagumo J, Arimoto S, Yoshizawa S (1964) An active pulse transmision line simulating nerve axon. Proc IRE 50:2061–2070 Ng YK, Lu X, Gulacsi A, Han W, Saxon MJ, Levitan ES (2003) Unexpected mobility variation among indyvidual secretory vesicles produces an apparent refractory neuropeptide pool. Biophys J 84:4127–4134 Oades RD, Daniels R, Rascher W (1998) Plasma neuropeptide—Y levels, monoamine metabolism, electro lyte excretion and drinking behavior in children with attention-deficit hyperactivity disorder. Psychiatry Res 80:177–186 Rall W (1969) Time constants and electronic length of membrane cylinders and neurons. Biophys J 9:1483–1508 Scott A (1995) Stairway to the mind. Springer, New York Scott R, Rusakov DA (2006) Main determinants of presynaptic Ca2+ dynamics at individual mossy fiber–CA3 pyramidal cell synapses. J Neurosci 26:7071–7081 Shakiryanova D, Tully A, Hewes RS, Deitcher DL, Levitan ES (2005) Activity-dependent liberation of synaptic neuropeptide vesicles. Nat Neurosci 8:173–178 Shakiryanova D, Tully A, Levitan ES (2006) Activity-dependent synaptic capture of transiting peptidergic vesicles. Nat Neurosci 9:896– 900 Tadeusiewicz R (1994) Problems of biocybernetics. PWN, Warszawa (in Polish) Virgo L, Humphries C, Mortimer A, Barnes T, Hirsch S, Belleroche Jde (1995) Cholecystokinin messenger RNA deficit in frontal and temporal cerebral cortex in schizophrenia. Biol Psychiatry 37: 694–701 Wong M, Shakiryanova D, Levitan ES (2008) Presynaptic ryanodine receptor-CamKII signaling is required for activity-dependent capture of transiting vesicles. J Mol Neurosci (to appear) doi:10.1007/ s12031-008-9080-8

123

458 Waxman SG (1972) Regional differentation of the axon: a review with special reference to the concept of the multiplex neuron. Brain Res 47:269–288

123

Biol Cybern (2008) 99:443–458 Texas Instruments Incorporated (2006) Precision ADC and DACs with 8051 microcontroller and flash memory (Rev. F). http://focus.ti. com/docs/prod/folders/print/msc1211y2.html, Texas Instruments